chemshell_qchem / blob
1 \r
2 #\r
3 # Transition-state search (with or without mode-following) or\r
4 # minimization with Baker algorithm\r
5 # "An algorithm for the location of transition states",\r
6 # J. Baker, J. Comp. Chem. 7 (1986) 385-395\r
7 #\r
8 # combination with DIIS\r
9 \r
10 proc newopt_baker2 opt {\r
11 \r
12 global no_parm\r
13 global no_log\r
14 global no_inh\r
15 \r
16 set collect_diis 0\r
17 set no_parm($opt,diis_dimension) 0\r
18 set diis_points 0\r
19 set operate_diis 0\r
20 set take_diis_step 0\r
21 if { $no_parm($opt,maxdiis) > [expr $no_parm($opt,nvar) - [llength $no_parm($opt,original_mode_list)] ] } {\r
22 set no_parm($opt,maxdiis) [expr $no_parm($opt,nvar) - [llength $no_parm($opt,original_mode_list)] ]\r
23 }\r
24 #\r
25 set restored_hessian 0\r
26 \r
27 set scale_step 1.0\r
28 \r
29 set no_parm($opt,nsteps) 0\r
30 # restart capability\r
31 switch $no_parm($opt,restart) {\r
32 1 {\r
33 puts stdout "newopt: attempt to restart Baker optimisation"\r
34 newopt_recover_restart_files $opt x f g h diis\r
35 $no lv\r
36 set no_log($opt,$no_parm($opt,nlog),f) [ get_matrix_element matrix=$opt.f indices= {0 0} format=%20.12f ]\r
37 set no_parm($opt,old_energy) $no_log($opt,$no_parm($opt,nlog),f)\r
38 #\r
39 if { $no_parm($opt,verbose) == 1} { print_matrix matrix=$opt.h }\r
40 #\r
41 }\r
42 0 {\r
43 #\r
44 # assumes a hessian is available\r
45 #\r
46 switch $no_parm($opt,input_hessian) {\r
47 undefined { chemerr "baker algorithm requires an input_hessian argument" }\r
48 default { }\r
49 }\r
50 $opt sh $no_parm($opt,input_hessian)\r
51 set no_parm($opt,hessian_ok) no\r
52 $opt fg\r
53 set no_parm($opt,old_energy) $no_log($opt,$no_parm($opt,nlog),f)\r
54 \r
55 if { $no_parm($opt,read_diis_space) } {\r
56 if { [ recover_diis_space ] } {\r
57 chemerr "could not read DIIS information"\r
58 }\r
59 }\r
60 }\r
61 }\r
62 \r
63 \r
64 \r
65 #\r
66 set no_parm($opt,new_energy) $no_parm($opt,old_energy)\r
67 puts stdout [format "$opt: initial energy: %14s" $no_parm($opt,old_energy) ]\r
68 set same_energy 0\r
69 \r
70 # save starting structure\r
71 $opt wr start\r
72 $opt pun\r
73 \r
74 if { $no_parm($opt,save_visualisation_files) } {\r
75 $opt wr vis\r
76 copy_object type=matrix from=$opt.h to=newopt.h_vis\r
77 }\r
78 \r
79 #\r
80 # set required hessian inertia \r
81 #\r
82 if { $no_parm($opt,follow_mode) == 0 } {\r
83 set no_parm($opt,pos) $no_parm($opt,nvar)\r
84 } else {\r
85 set no_parm($opt,pos) [expr $no_parm($opt,nvar) - 1 ]\r
86 }\r
87 # set no_parm($opt,best_hessian) $no_parm($opt,nvar)\r
88 set no_parm($opt,best_step) 3\r
89 #\r
90 # initialize the loop\r
91 #\r
92 set forced_conv 0\r
93 set more 1\r
94 # diagonalise hessian\r
95 $opt dh\r
96 # check hessian character\r
97 $opt hc\r
98 set no_inh($opt,pos_before) $no_inh($opt,pos)\r
99 set no_parm($opt,best_hessian) [expr $no_parm($opt,pos) - $no_inh($opt,pos) ]\r
100 set no_parm($opt,hessian_update_old) $no_parm($opt,hessian_update)\r
101 #\r
102 # set original hessian eigenmodes\r
103 #\r
104 switch $no_parm($opt,original_modes) {\r
105 undefined { set no_parm($opt,original_modes) $opt.evech }\r
106 default { } \r
107 }\r
108 puts stdout [format "newopt: original_modes=%s" $no_parm($opt,original_modes) ]\r
109 $opt sms $no_parm($opt,original_modes)\r
110 #\r
111 # start search\r
112 #\r
113 while {$more} {\r
114 \r
115 #\r
116 # if previous hessian is used, reset position and gradient to\r
117 # that belonging to the hessian for hessian updates\r
118 #\r
119 if { $restored_hessian } {\r
120 copy_object type=matrix from=$opt.px to=$opt.x\r
121 copy_object type=matrix from=$opt.pg to=$opt.g\r
122 set restored_hessian 0\r
123 }\r
124 #\r
125 set go_back 0\r
126 if {$no_parm($opt,follow_mode) != 0} {\r
127 #\r
128 # this is a TS search; set or select the required mode\r
129 #\r
130 if { $no_parm($opt,freeze_modes) } {\r
131 catch { unset no_parm($opt,map_modes) }\r
132 set sms_result [$opt sms]\r
133 switch $sms_result {\r
134 0 {\r
135 set new_mode_list ""\r
136 foreach el $no_parm($opt,original_mode_list) {\r
137 lappend new_mode_list [lindex $no_parm($opt,map_modes) [expr $el - 1]]\r
138 }\r
139 }\r
140 default {\r
141 set fmode 1\r
142 set mmode 1\r
143 foreach el $no_parm($opt,original_mode_list) {\r
144 if { $el == $no_parm($opt,follow_mode) } { set fmode 0 }\r
145 if { $el == $sms_result } { set mmode 0 }\r
146 }\r
147 if { $fmode } {\r
148 puts stdout "newopt: strong mixing in TS mode: go back and freeze TS mode"\r
149 lappend no_parm($opt,original_mode_list) $no_parm($opt,follow_mode)\r
150 set go_back 1\r
151 }\r
152 if { $mmode } {\r
153 puts stdout [format "newopt: strong mixing in TS mode: go back and freeze mode %s" $sms_result ]\r
154 lappend no_parm($opt,original_mode_list) $sms_result\r
155 set go_back 1\r
156 }\r
157 set new_mode_list $no_parm($opt,original_mode_list)\r
158 if { $go_back } {\r
159 copy_object type=matrix from=$opt.px to=$opt.x\r
160 copy_object type=matrix from=$opt.pg to=$opt.g\r
161 copy_object type=matrix from=$opt.ph to=$opt.h\r
162 $opt dh\r
163 $opt hc\r
164 }\r
165 }\r
166 }\r
167 set no_parm($opt,mode_list) $new_mode_list\r
168 }\r
169 \r
170 #\r
171 $opt sm\r
172 #\r
173 # refresh the TS mode-following vector if required\r
174 #\r
175 if {$no_parm($opt,update_mode)} {$opt um}\r
176 \r
177 if {([expr $no_inh($opt,neg) + $no_inh($opt,zero)] != 1)} {\r
178 # do not count zero-value eigenmodes if they are to be ignored\r
179 if { $no_parm($opt,freeze_soft_modes) } {\r
180 #\r
181 if {$no_inh($opt,neg) != 1} {\r
182 set no_parm($opt,hessian_ok) no\r
183 puts stdout "the hessian does not have the required inertia"\r
184 }\r
185 } else {\r
186 # the hessian does not have the required inertia\r
187 # OR the selected mode is not the lowest\r
188 #\r
189 set no_parm($opt,hessian_ok) no\r
190 puts stdout "the hessian does not have the required inertia"\r
191 }\r
192 if { $no_parm($opt,condition_hessian) } {\r
193 #no_gn\r
194 #if { $no_parm($opt,gnorm) > $no_parm($opt,conv_rmsgrad) } \r
195 puts stdout "hessian conditioning"\r
196 $opt mh\r
197 $opt hc\r
198 }\r
199 } elseif {$no_parm($opt,follow_mode) != 1} {\r
200 \r
201 puts stdout "no_parm($opt,follow_mode) = $no_parm($opt,follow_mode)"\r
202 \r
203 puts stdout "the hessian has the correct inertia, but not the correct structure"\r
204 #\r
205 # just to try and make it more robust\r
206 # follow the lowest mode if it is the only negative eigenvalue one\r
207 # and use it for the rest of the search\r
208 #\r
209 # set no_parm($opt,hessian_ok) yes\r
210 set no_parm($opt,hessian_ok) no\r
211 if { $no_parm($opt,condition_hessian) } {\r
212 puts stdout "hessian conditioning"\r
213 $opt mh\r
214 $opt hc\r
215 }\r
216 #\r
217 # the fact that the mode to follow isn't the lowest one \r
218 # should not deter convergence; that is why\r
219 # the hessian is flagged ok\r
220 # be sure to check, though - and start again if\r
221 # not satisfied with the result\r
222 #\r
223 # set no_parm($opt,follow_mode) 1\r
224 # no_um\r
225 } else {\r
226 set no_parm($opt,hessian_ok) yes\r
227 puts stdout "the hessian has the correct structure"\r
228 }\r
229 } else {\r
230 #\r
231 # minimization\r
232 #\r
233 if {[expr $no_inh($opt,pos) + $no_inh($opt,zero)] == $no_parm($opt,nvar)} {\r
234 if { $no_parm($opt,freeze_soft_modes) || ($no_inh($opt,zero) == 0) } {\r
235 set no_parm($opt,hessian_ok) yes\r
236 puts stdout "the hessian has the correct structure"\r
237 } else {\r
238 set no_parm($opt,hessian_ok) no\r
239 puts stdout "the hessian does not have the correct structure"\r
240 if {$no_parm($opt,condition_hessian) } {\r
241 puts stdout "hessian conditioning"\r
242 $opt mh\r
243 $opt hc\r
244 }\r
245 }\r
246 } else {\r
247 set no_parm($opt,hessian_ok) no\r
248 if {$no_parm($opt,condition_hessian) } {\r
249 # no_gn\r
250 # if { $no_parm($opt,gnorm) > $no_parm($opt,conv_rmsgrad) } \r
251 puts stdout "hessian conditioning"\r
252 $opt mh\r
253 $opt hc\r
254 }\r
255 }\r
256 }\r
257 \r
258 #\r
259 if { $no_parm($opt,nsteps) > 2 } {\r
260 #\r
261 # calculate overlap with previous step\r
262 #\r
263 transpose_matrix matrix=$opt.ps result=$opt.pst\r
264 multiply_matrix left=$opt.pst right=$opt.s result=$opt.olaps\r
265 multiply_matrix left=$opt.ts right=$opt.s result=$opt.olapss\r
266 multiply_matrix left=$opt.pst right=$opt.ps result=$opt.olapp\r
267 set ovlap [get_matrix_element matrix=$opt.olaps indices= {0 0} format=%e]\r
268 set ovlass [get_matrix_element matrix=$opt.olapss indices= {0 0} format=%e]\r
269 set ovlapp [get_matrix_element matrix=$opt.olapp indices= {0 0} format=%e]\r
270 \r
271 if { ($ovlass > 1.e-10) && ($ovlapp > 1.e-10) } {\r
272 puts stdout [format "newopt: overlap with previous step %f" [expr $ovlap/(sqrt([expr $ovlass*$ovlapp]))] ]\r
273 }\r
274 if { ($no_parm($opt,hessian_ok) == "yes") && (!$go_back) } {\r
275 #\r
276 # collect points for diis extrapolation in correct hessian region\r
277 #\r
278 set collect_diis 1\r
279 } else {\r
280 set collect_diis 0\r
281 }\r
282 }\r
283 \r
284 #\r
285 # PAUL - need to integrate the linear dependency test\r
286 # and the decision to remove points\r
287 # see block - $diis_points < $no_parm($opt,mindiis)\r
288 \r
289 \r
290 if { $collect_diis } { \r
291 if { $no_parm($opt,diis_dimension) >= $no_parm($opt,maxdiis) } {\r
292 \r
293 if { $diis_points < $no_parm($opt,mindiis) } {\r
294 #\r
295 # remove earliest point\r
296 #\r
297 $opt gdd 0\r
298 } else {\r
299 #\r
300 # remove point with largest RMS gradient\r
301 #\r
302 $opt gddg\r
303 }\r
304 incr no_parm($opt,diis_dimension) -1\r
305 }\r
306 if { [ $opt gda ] == 0 } then {\r
307 if { $no_parm($opt,diis_print) == 1 } { $opt gdp }\r
308 incr no_parm($opt,diis_dimension) \r
309 incr diis_points\r
310 }\r
311 \r
312 if { ($diis_points >= $no_parm($opt,diis_steps)) && \r
313 ($no_parm($opt,diis_dimension) >= $no_parm($opt,mindiis)) && \r
314 $no_parm($opt,diis_baker)} {\r
315 set operate_diis 1\r
316 } \r
317 }\r
318 \r
319 #\r
320 # before doing so, check if any modes have changed order and set options accordingly\r
321 #\r
322 if { $no_parm($opt,go_with_diis) && \r
323 ($no_parm($opt,diis_dimension) >= $no_parm($opt,mindiis)) && \r
324 ($no_parm($opt,hessian_ok) == "yes") } {\r
325 \r
326 copy_object type=matrix from=$opt.x to=$opt.x.tmp\r
327 copy_object type=matrix from=$opt.g to=$opt.g.tmp\r
328 #\r
329 # invert hessian\r
330 #\r
331 $opt ih\r
332 #\r
333 # calculate diis step\r
334 #\r
335 $opt gds\r
336 #\r
337 if { $no_parm($opt,verbose) == 1 } {\r
338 puts stdout "====================== diis extrapolations ====================="\r
339 print_matrix matrix=$opt.g\r
340 $opt grms\r
341 puts stdout [format "newopt: predicted GDIIS gradient rms: %f" $no_parm($opt,grms)]\r
342 print_matrix matrix=$opt.x\r
343 puts stdout "================================================================"\r
344 }\r
345 #\r
346 $opt fg\r
347 $opt grms\r
348 puts stdout [ format "newopt: gdiis gradient rms %f" $no_parm($opt,grms) ]\r
349 set no_parm($opt,new_energy) $no_log($opt,$no_parm($opt,nlog),f)\r
350 set diff [expr $no_parm($opt,new_energy) - $no_parm($opt,old_energy) ]\r
351 puts stdout [format "$opt: new energy: %14s, energy change: %14s" $no_log($opt,$no_parm($opt,nlog),f) $diff ]\r
352 set no_parm($opt,old_energy) $no_parm($opt,new_energy)\r
353 #\r
354 subtract_matrix left=$opt.x right=$opt.x.tmp result=$opt.s\r
355 copy_object type=matrix from=$opt.x.tmp to=$opt.px\r
356 copy_object type=matrix from=$opt.g.tmp to=$opt.pg\r
357 switch $no_parm($opt,hessian_update) {\r
358 bfgs { $opt bfu }\r
359 powell { $opt pu }\r
360 mursar { $opt msu }\r
361 hybrid { $opt pmsu }\r
362 fde { $opt efdh $no_parm($opt,var_list) }\r
363 fde2 { $opt efdh2 $no_parm($opt,var_list) }\r
364 fd { $opt fdh }\r
365 fd2 { $opt fdh2 }\r
366 none { }\r
367 default { \r
368 chemerr "Baker: illegal hessian_update argument $no_parm($opt,hessian_update)"\r
369 }\r
370 }\r
371 $opt dh\r
372 $opt hc\r
373 \r
374 #\r
375 # PAUL - how does no_parm($opt,hessian_ok) get set\r
376 # this test is out of date\r
377 #\r
378 \r
379 \r
380 if { $no_parm($opt,hessian_ok) == "yes" } {\r
381 if { $no_parm($opt,diis_dimension) >= $no_parm($opt,maxdiis) } {\r
382 #\r
383 # PAUL - need to check on meaning of the following test\r
384 #\r
385 if { $diis_points < $no_parm($opt,mindiis) } {\r
386 #\r
387 # remove earliest point\r
388 #\r
389 $opt gdd 0\r
390 } else {\r
391 #\r
392 # remove point with largest RMS gradient\r
393 #\r
394 $opt gddg\r
395 }\r
396 incr no_parm($opt,diis_dimension) -1\r
397 }\r
398 if { [ $opt gda ] == 0 } then {\r
399 if { $no_parm($opt,diis_print) == 1 } { $opt gdp }\r
400 incr no_parm($opt,diis_dimension)\r
401 }\r
402 } \r
403 }\r
404 \r
405 #\r
406 # --- end of diis section (go_with_diis alternating prfo/diis strategy) -----\r
407 #\r
408 \r
409 #\r
410 # Compute PRFO step\r
411 #\r
412 # transform the gradient matrix with transpose of eigenvector matrix\r
413 $opt ug\r
414 \r
415 # calculate the step (with mode-following if required)\r
416 switch [ $opt prfo ] {\r
417 0 {\r
418 #\r
419 # reverse step-vector if follow_mode was negative\r
420 #\r
421 if {$no_parm($opt,follow_mode) < 0} {\r
422 upscale_matrix matrix=$opt.s factor=-1.0 result=$opt.s \r
423 }\r
424 $opt cl $no_parm($opt,maxlength)\r
425 #\r
426 # an extra step scale factor may be used to\r
427 # deal with poor hessian and failed update\r
428 #\r
429 if { $scale_step != 1.0 } {\r
430 upscale_matrix matrix=$opt.s result=$opt.s factor=$scale_step\r
431 set scale_step 1.0\r
432 }\r
433 #\r
434 $opt s\r
435 #\r
436 if { $no_parm($opt,verbose) == 1 } { print_matrix matrix=$opt.s }\r
437 if {$no_parm($opt,list_variables) } {$opt lv}\r
438 #\r
439 # estimate the energy change in the calculated step\r
440 #\r
441 $opt ee;\r
442 #\r
443 # abort the job if maximum number of steps is set to 0\r
444 #\r
445 if { $no_parm($opt,maxstep) == 0 } {\r
446 set more 0\r
447 #\r
448 # not sure what to do here\r
449 #set no_parm($opt,status) -1\r
450 } else {\r
451 #\r
452 # calculate the function value at the new point\r
453 #\r
454 $opt f\r
455 #\r
456 # print progress of the optimization\r
457 #\r
458 #\r
459 # PAUL - I don't think we should allow the algorithm to be sensitive\r
460 # to the log data?\r
461 # use get_matrix_element matrix=f ..\r
462 # Why should we need to repeat the energy calc\r
463 #\r
464 set no_parm($opt,new_energy) $no_log($opt,$no_parm($opt,nlog),f)\r
465 if {$no_parm($opt,new_energy) == " "} {\r
466 $opt f\r
467 set no_parm($opt,new_energy) $no_log($opt,$no_parm($opt,nlog),f)\r
468 if {$no_parm($opt,new_energy) == " "} {\r
469 chemerr "something wrong with energy calculation"\r
470 }\r
471 }\r
472 set diff [expr $no_parm($opt,new_energy) - $no_parm($opt,old_energy) ]\r
473 puts stdout [format "$opt: new energy: %14s, energy change: %14s" $no_log($opt,$no_parm($opt,nlog),f) $diff ]\r
474 if {($diff < $no_parm($opt,toler_energy)) && ($diff > -$no_parm($opt,toler_energy)) } {\r
475 #\r
476 # The energy has not changed within the given threshold\r
477 # When this happens more than once, abort the search\r
478 #\r
479 incr same_energy 1\r
480 } else {\r
481 set same_energy 0\r
482 }\r
483 if { ($same_energy > 0) } {\r
484 set forced_conv 1\r
485 puts stdout [format "newopt: forced convergence because energy change is smaller than threshold: %e twice in a row" $no_parm($opt,toler_energy)] \r
486 }\r
487 if { $forced_conv } {\r
488 #\r
489 # One exception: When modes are frozen, now free them up and continue the search \r
490 #\r
491 if { $no_parm($opt,freeze_modes) && $no_parm($opt,mode_list) != "" && $no_parm($opt,release_modes) != 0 } {\r
492 puts stdout "newopt: continuing search with less frozen modes"\r
493 set frozen [ llength $no_parm($opt,mode_list) ]\r
494 set eigenmode_list ""\r
495 for { set i 0 } { $i < [expr $frozen - $no_parm($opt,release_modes)] } { incr i 1 } {\r
496 lappend eigenmode_list [ lindex $no_parm($opt,mode_list) $i ]\r
497 }\r
498 set no_parm($opt,mode_list) $eigenmode_list\r
499 set frozen [ llength $no_parm($opt,original_mode_list) ]\r
500 set eigenmode_list ""\r
501 for { set i 0 } { $i < [expr $frozen - $no_parm($opt,release_modes)] } { incr i 1 } {\r
502 puts stdout [format "newopt: mode %d frozen" [expr $i + 1]]\r
503 lappend eigenmode_list [ lindex $no_parm($opt,original_mode_list) $i ]\r
504 }\r
505 set no_parm($opt,original_mode_list) $eigenmode_list\r
506 if { $no_parm($opt,update_mode) } { copy_object type=matrix from=evech to=modes }\r
507 set same_energy 0\r
508 set no_parm($opt,hessian_update) $no_parm($opt,hessian_update_old)\r
509 \r
510 set forced_conv 0\r
511 } else {\r
512 set more 0\r
513 }\r
514 }\r
515 #\r
516 # with small energy changes, do not update the hessian any further\r
517 #\r
518 if {($diff < $no_parm($opt,toler_no_update)) && \\r
519 ($diff > -$no_parm($opt,toler_no_update)) && \\r
520 ($no_parm($opt,hessian_ok) == "yes")} {\r
521 \r
522 puts stdout "newopt: hessian no longer updated due to small energy difference from last step "\r
523 set no_parm($opt,hessian_update) none\r
524 }\r
525 #\r
526 # recalculate maximum step length based on experience\r
527 #\r
528 switch $no_parm($opt,trust_region) {\r
529 active {\r
530 set more2 1\r
531 set npass 0 \r
532 copy_object type=matrix from=$opt.x to=$opt.x.tmp\r
533 copy_object type=matrix from=$opt.g to=$opt.g.tmp\r
534 while {$more2} {\r
535 incr npass\r
536 switch [ $opt tr ] {\r
537 1 { \r
538 #\r
539 # hessian predicted the wrong sign,\r
540 # or the ratio between predicted and actual energy difference is very high or low:\r
541 # rescale the step and step again\r
542 #\r
543 # fi { $no_parm($opt,maxlength) > $no_parm($opt,stepmin) } \r
544 # in a minimization, if the energy is lower than predicted, take the step anyways \r
545 \r
546 # if { ($follow_mode == 0) && ( $diff < [get_matrix_element matrix=dele indices= {0 0} format=%e]) } { set more2 0 } else { close brace added for indent -> }\r
547 \r
548 \r
549 $opt cl $no_parm($opt,maxlength)\r
550 $opt s\r
551 if { $no_parm($opt,verbose) == 1 } { print_matrix matrix=$opt.s }\r
552 $opt ee\r
553 $opt f\r
554 set no_parm($opt,new_energy) $no_log($opt,$no_parm($opt,nlog),f)\r
555 set diff [expr $no_parm($opt,new_energy) - $no_parm($opt,old_energy) ]\r
556 puts stdout [format "$opt: new energy: %14s, energy change: %14s" $no_log($opt,$no_parm($opt,nlog),f) $diff ]\r
557 # close brace was here\r
558 if { $npass > 1 } { set more2 0 }\r
559 }\r
560 0 {\r
561 set more2 0 \r
562 }\r
563 default { \r
564 set more2 0 \r
565 }\r
566 } \r
567 if {($npass > 4) || ($no_parm($opt,maxlength) == $no_parm($opt,stepmin))} {\r
568 set more2 0\r
569 } \r
570 }\r
571 if { $npass > 1 } {\r
572 copy_object type=matrix from=$opt.x.tmp to=$opt.px\r
573 copy_object type=matrix from=$opt.g.tmp to=$opt.pg\r
574 } \r
575 }\r
576 default { }\r
577 }\r
578 $opt pun\r
579 \r
580 if { $no_parm($opt,save_visualisation_files) } {\r
581 $opt wr vis\r
582 copy_object type=matrix from=$opt.h to=newopt.h_vis\r
583 }\r
584 \r
585 #\r
586 # PAUL - question of diis vs prfo\r
587 #\r
588 # calculate the gradient at the new point\r
589 $opt g\r
590 #\r
591 if { $operate_diis } {\r
592 #\r
593 $opt grms\r
594 set prfo_grms $no_parm($opt,grms)\r
595 puts stdout [ format "newopt: prfo gradient rms %f" $prfo_grms ]\r
596 #\r
597 #\r
598 # PAUL - add save of f here\r
599 # we need to document the use of f/pf etc\r
600 #\r
601 # save the current and previous position/function/gradient\r
602 #\r
603 copy_object type=matrix from=$opt.x to=$opt.x.tmp\r
604 copy_object type=matrix from=$opt.g to=$opt.g.tmp\r
605 copy_object type=matrix from=$opt.f to=$opt.f.tmp\r
606 copy_object type=matrix from=$opt.px to=$opt.px.tmp\r
607 copy_object type=matrix from=$opt.pg to=$opt.pg.tmp\r
608 #\r
609 # invert hessian\r
610 #\r
611 $opt ih\r
612 #\r
613 # calculate diis step\r
614 #\r
615 #\r
616 # PAUL - failure to solve diis eqns is fatal?\r
617 # maybe shouldn't be\r
618 $opt gds\r
619 #\r
620 if { $no_parm($opt,verbose) == 1 } {\r
621 puts stdout "====================== diis extrapolations ====================="\r
622 print_matrix matrix=$opt.g\r
623 print_matrix matrix=$opt.x\r
624 puts stdout "================================================================"\r
625 }\r
626 #\r
627 \r
628 #\r
629 # PAUL - maybe we can give up on diis if extrapolation fails\r
630 #\r
631 $opt grms\r
632 puts stdout [ format "newopt: extrapolated gdiis gradient rms %f" $no_parm($opt,grms) ]\r
633 \r
634 set take_diis_step 0\r
635 #\r
636 # PAUL - apply test on extrapolated RMS\r
637 #\r
638 if { $no_parm($opt,grms) < $prfo_grms } {\r
639 $opt fg\r
640 $opt grms\r
641 puts stdout [ format "newopt: gdiis gradient rms %f" $no_parm($opt,grms) ]\r
642 #\r
643 # test gradient norm of DIIS step against prfo step\r
644 #\r
645 if { $no_parm($opt,grms) > $prfo_grms } {\r
646 # go with prfo step\r
647 puts stdout "newopt: prfo step retained: prfo grms $prfo_grms is better than diis rms $no_parm($opt,grms)"\r
648 } else {\r
649 puts stdout "newopt: GDIIS retained: prfo grms $prfo_grms is poorer than diis rms $no_parm($opt,grms)"\r
650 subtract_matrix left=$opt.x right=$opt.px.tmp result=$opt.s\r
651 set take_diis_step 1\r
652 copy_object type=matrix from=$opt.px.tmp to=$opt.px\r
653 copy_object type=matrix from=$opt.pg.tmp to=$opt.pg\r
654 }\r
655 } else {\r
656 puts stdout "newopt: prfo step retained: prfo grms $prfo_grms is better than extrapolated diis rms $no_parm($opt,grms)"\r
657 }\r
658 \r
659 if { ! $take_diis_step } {\r
660 copy_object type=matrix from=$opt.f.tmp to=$opt.f\r
661 copy_object type=matrix from=$opt.x.tmp to=$opt.x\r
662 copy_object type=matrix from=$opt.g.tmp to=$opt.g\r
663 }\r
664 }\r
665 # test maximum number of steps\r
666 \r
667 #\r
668 # PAUL move this up to avoid unnecessary hessian evaluations\r
669 #\r
670 set converged [ $opt tc ]\r
671 \r
672 if { $no_parm($opt,nsteps) >= $no_parm($opt,maxstep) } {\r
673 puts stdout [ format "maximum number of steps %d exceeded" $no_parm($opt,maxstep) ]\r
674 $opt wr best\r
675 set more 0\r
676 set no_parm($opt,status) -1\r
677 } else {\r
678 #\r
679 # reset diis controls\r
680 #\r
681 if { $operate_diis } {\r
682 set operate_diis 0\r
683 #\r
684 # if diis step was accepted give diis another chance\r
685 # otherwise start counting again for next try\r
686 #\r
687 if { !$take_diis_step } { set diis_points 0 }\r
688 }\r
689 \r
690 #\r
691 \r
692 set do_update [expr fmod($no_parm($opt,nsteps),$no_parm($opt,recompute_hessian))]\r
693 \r
694 puts stdout "do update $do_update converged $converged $no_parm($opt,hessian_update)"\r
695 \r
696 # we may have convergence, but with an incorrect hessian!\r
697 # therefore, we must try the hessian update\r
698 # fi open ! $converged close open\r
699 \r
700 if { $do_update } {\r
701 \r
702 set errcode 0\r
703 #\r
704 # PAUL - are some of these sensible strategies - fde and fde2 only modify some of the \r
705 # hessian, perhaps they should be combined with true updates?\r
706 #\r
707 # nb - need to make special provisions for diis solutions\r
708 #\r
709 #\r
710 switch $no_parm($opt,hessian_update) {\r
711 bfgs { $opt bfu }\r
712 powell { set errcode [$opt pu] }\r
713 mursar { set errcode [ $opt msu ] }\r
714 hybrid { set errcode [$opt pmsu] }\r
715 fde { $opt efdh $no_parm($opt,var_list) }\r
716 fde2 { $opt efdh2 $no_parm($opt,var_list) }\r
717 fd { $opt fdh }\r
718 fd2 { $opt fdh2 }\r
719 none { }\r
720 default { chemerr "Baker: illegal hessian_update argument $no_parm($opt,hessian_update)" }\r
721 }\r
722 \r
723 if { $errcode != 0 } {\r
724 \r
725 \r
726 switch $no_parm($opt,update_failure_strategy) {\r
727 \r
728 ignore {\r
729 puts stdout "the old hessian will be retained"\r
730 }\r
731 \r
732 efdh_max {\r
733 if { $no_parm($opt,follow_mode) } {\r
734 set fdelist "maxgrad maxstep $no_parm($opt,follow_mode)"\r
735 } else {\r
736 set fdelist "maxgrad maxstep"\r
737 }\r
738 $opt efdh $fdelist\r
739 }\r
740 \r
741 upscale_step {\r
742 puts stdout "next step will be scaled by 3.0"\r
743 set scale_step 3.0\r
744 }\r
745 \r
746 recompute_hessian {\r
747 puts stdout "Full hessian recomputation forced by failed update"\r
748 $opt h \r
749 }\r
750 \r
751 default {\r
752 chemerr "newopt: argument update_failure_strategy incorrectly set"\r
753 }\r
754 }\r
755 }\r
756 } else {\r
757 #\r
758 # hessian recalculation\r
759 #\r
760 switch $no_parm($opt,recompute_method) {\r
761 fd { $opt fdh }\r
762 fd2 { $opt fdh2 }\r
763 analytic { $opt h }\r
764 fde { $opt efdh $no_parm($opt,var_list) }\r
765 fde2 { $opt efdh2 $no_parm($opt,var_list) }\r
766 none {\r
767 copy_object type=matrix from=$opt.ph to=$opt.h\r
768 set restored_hessian 1\r
769 }\r
770 ignore { }\r
771 default {chemerr " Baker: illegal recompute_method argument " }\r
772 }\r
773 }\r
774 \r
775 # diagonalise hessian\r
776 $opt dh\r
777 # check hessian character\r
778 $opt hc\r
779 \r
780 # close\r
781 # close\r
782 }\r
783 \r
784 # test convergence\r
785 \r
786 if { $no_parm($opt,freeze_modes) } {\r
787 #\r
788 # project out contributions of frozen modes to the gradient\r
789 # \r
790 copy_object type=matrix from=$opt.g to=$opt.g.tmp\r
791 $opt eg \r
792 }\r
793 \r
794 if { ( $converged == 1) && ($no_parm($opt,hessian_ok) == "yes") } { \r
795 \r
796 if { $no_parm($opt,freeze_modes) && $no_parm($opt,mode_list) != "" && $no_parm($opt,release_modes) != 0 } {\r
797 #\r
798 # Normally this is convergence, but in case of frozen modes\r
799 # we need to continue after freeing them up\r
800 #\r
801 # Here, the indices of the frozen modes are reset and the list shortened\r
802 #\r
803 set frozen [ llength $no_parm($opt,mode_list) ]\r
804 set eigenmode_list ""\r
805 for { set i 0 } { $i < [expr $frozen - $no_parm($opt,release_modes)] } { incr i 1 } {\r
806 lappend eigenmode_list [ lindex $no_parm($opt,mode_list) $i ]\r
807 }\r
808 set no_parm($opt,mode_list) $eigenmode_list\r
809 set frozen [ llength $no_parm($opt,original_mode_list) ]\r
810 set eigenmode_list ""\r
811 for { set i 0 } { $i < [expr $frozen - $no_parm($opt,release_modes)] } { incr i 1 } {\r
812 lappend eigenmode_list [ lindex $no_parm($opt,original_mode_list) $i ]\r
813 }\r
814 set no_parm($opt,original_mode_list) $eigenmode_list\r
815 if { $no_parm($opt,update_mode) } { copy_object type=matrix from=evech to=modes }\r
816 \r
817 incr no_parm($opt,mindiis) $no_parm($opt,release_modes)\r
818 incr no_parm($opt,maxdiis) $no_parm($opt,release_modes)\r
819 if { $no_parm($opt,maxdiis) > [expr $no_parm($opt,nvar) - [llength $no_parm($opt,original_mode_list)] ] } {\r
820 set no_parm($opt,maxdiis) [expr $no_parm($opt,nvar) - [llength $no_parm($opt,original_mode_list)] ]\r
821 }\r
822 } else {\r
823 set more 0\r
824 }\r
825 }\r
826 # restore the full gradient\r
827 if { $no_parm($opt,freeze_modes) } { \r
828 copy_object type=matrix from=$opt.g.tmp to=$opt.g\r
829 }\r
830 # see if the current point makes sense\r
831 if {$no_parm($opt,run_diagnostics)} {\r
832 if { [$opt d] } {\r
833 # abort the search - there's something wrong\r
834 set more 0\r
835 set no_parm($opt,status) -1\r
836 } else {\r
837 #\r
838 # check for changes in the hessian inertia\r
839 # if these occur, undo the updating of the hessian\r
840 # unless it was reached through finite differencing along \r
841 # eigenmodes or recomputation of parts of the hessian\r
842 #\r
843 if { [expr $no_inh($opt,pos_before) - $no_inh($opt,pos) ] >= 1 } {\r
844 #\r
845 # check if this is a TS search and if we have one negative eigenvalue\r
846 #\r
847 if { $no_parm($opt,follow_mode) != 0 } {\r
848 if { !($no_inh($opt,neg) == 1) } {\r
849 #\r
850 # we do not have the correct number of negative eigenvalues\r
851 # first try other update\r
852 #\r
853 if { ($no_parm($opt,hessian_update) != "fde") && [expr fmod($no_parm($opt,nsteps),$no_parm($opt,recompute_hessian))] } {\r
854 set no_inh($opt,pos_other_update) $no_inh($opt,pos)\r
855 copy_object type=matrix from=$opt.h to=$opt.h_other\r
856 copy_object type=matrix from=$opt.ph to=$opt.h\r
857 #\r
858 # PAUL\r
859 # hybrid { $opt pu }\r
860 switch $no_parm($opt,hessian_update) {\r
861 mursar { $opt pmsu }\r
862 powell { $opt pmsu }\r
863 hybrid {\r
864 puts stdout "apparent hessian degradation detected - recompute"\r
865 $opt efdh $no_parm($opt,var_list)\r
866 }\r
867 default { }\r
868 }\r
869 $opt dh\r
870 $opt hc\r
871 # check for improvement\r
872 if { [expr $no_inh($opt,pos_before) - $no_inh($opt,pos) ] >= 1 } {\r
873 # we've still got a worse hessian than before\r
874 if { $no_parm($opt,fin_diff_eigenmodes) } {\r
875 # finite-difference selected eigenmodes\r
876 set eigenmode_list ""\r
877 for { set i 1 } { $i <= [expr $no_inh($opt,neg) + $no_inh($opt,zero)] } { incr i 1 } {\r
878 lappend eigenmode_list $i\r
879 }\r
880 puts stdout "newopt: finite differencing of selected eigenmodes to improve the hessian"\r
881 $opt efdh2 $eigenmode_list\r
882 $opt dh\r
883 $opt hc\r
884 }\r
885 }\r
886 }\r
887 # check for improvement\r
888 if { [expr $no_inh($opt,pos_before) - $no_inh($opt,pos) ] >= 1 } {\r
889 # we've still got a worse hessian than before\r
890 if { $no_parm($opt,keep_correct_hessian) } {\r
891 # restore previous hessian\r
892 puts stdout "newopt: hessian updates ignored --- previous hessian restored"\r
893 copy_object type=matrix from=$opt.ph to=$opt.h\r
894 set restored_hessian 1\r
895 $opt dh\r
896 $opt hc\r
897 } else {\r
898 # revert to other update if that was better\r
899 if { $no_inh($opt,pos) < $no_inh($opt,pos_other_update) } {\r
900 puts stdout "newopt: first hessian update is better, therefore restored"\r
901 copy_object type=matrix from=$opt.h_other to=$opt.h\r
902 $opt dh\r
903 $opt hc\r
904 }\r
905 }\r
906 }\r
907 }\r
908 } \r
909 } else {\r
910 if { ($no_parm($opt,hessian_ok) == "yes") && ( [expr $no_inh($opt,pos_before) - $no_inh($opt,pos) ] == -1) } {\r
911 #\r
912 # hessian was ok but has now gone off the required inertia,\r
913 # showing positive eigenvalues only \r
914 # restore previous hessian\r
915 #\r
916 puts stdout "newopt: hessian updates ignored --- previous hessian restored"\r
917 copy_object type=matrix from=$opt.ph to=$opt.h\r
918 set restored_hessian 1\r
919 $opt dh\r
920 $opt hc\r
921 }\r
922 }\r
923 # save the point for a restart if better suited than any of the previous points\r
924 # presently, the best restart structure is determined on the basis of the lowest rms gradient. \r
925 # only start saving after 3 steps have been taken, in case the search started at a stationary point\r
926 if { $no_parm($opt,nsteps) > 3 } { \r
927 # check for best restart point \r
928 # making sure the hessian inertia is the best approximation to the type of point we're after\r
929 if { ($no_inh($opt,pos) <= $no_parm($opt,pos)) && ( [expr $no_parm($opt,pos) - $no_inh($opt,pos) ] <= $no_parm($opt,best_hessian)) } {\r
930 if { [$opt sp] } {\r
931 set no_parm($opt,best_hessian) [expr $no_parm($opt,pos) - $no_inh($opt,pos) ]\r
932 set no_parm($opt,best_step) $no_parm($opt,nsteps)\r
933 if {$no_parm($opt,set_restart)} {\r
934 copy_object type=matrix from=$opt.h to=newopt.h_best\r
935 $opt wr best\r
936 }\r
937 }\r
938 }\r
939 puts stdout "check steps $no_parm($opt,best_step) $nsteps $no_parm($opt,trust_steps) "\r
940 if {( [expr $no_parm($opt,nsteps) - $no_parm($opt,best_step) ] > $no_parm($opt,trust_steps)) } {\r
941 # no improvement found in the last $no_parm($opt,trust_steps) steps; abort the search for a restart\r
942 puts stdout [ format "newopt: search abandoned because last %d steps did not bring improvement" $no_parm($opt,trust_steps) ]\r
943 $opt sr\r
944 set more 0\r
945 set no_parm($opt,status) -1\r
946 }\r
947 }\r
948 }\r
949 }\r
950 incr no_parm($opt,nsteps)\r
951 }\r
952 }\r
953 1 {\r
954 # the p-rfo step has failed\r
955 if { $no_parm($opt,diis_baker) && ($no_parm($opt,diis_dimension) >= $no_parm($opt,mindiis)) } {\r
956 puts stdout "newopt: perform diis extrapolation after p-rfo failure "\r
957 #\r
958 # invert hessian\r
959 #\r
960 $opt ih\r
961 #\r
962 # calculate diis step\r
963 #\r
964 $opt gds\r
965 #\r
966 if { $no_parm($opt,verbose) == 1 } {\r
967 puts stdout "====================== diis extrapolations ====================="\r
968 print_matrix matrix=$opt.g\r
969 print_matrix matrix=$opt.x\r
970 puts stdout "================================================================"\r
971 }\r
972 #\r
973 subtract_matrix left=$opt.x right=$opt.px result=$opt.s\r
974 $opt fg\r
975 switch $no_parm($opt,hessian_update) {\r
976 bfgs { $opt bfu }\r
977 powell { $opt pu }\r
978 mursar { $opt msu }\r
979 hybrid { $opt pmsu }\r
980 fde { $opt efdh $no_parm($opt,var_list) }\r
981 fde2 { $opt efdh2 $no_parm($opt,var_list) }\r
982 none { }\r
983 default { chemerr "Baker: illegal hessian_update argument $no_parm($opt,hessian_update)" }\r
984 }\r
985 } else {\r
986 #\r
987 # the p-rfo failure may be due to a bad hessian\r
988 # recalculate the hessian, or take the previous one,\r
989 # and try the step again, reducing the \r
990 # size of the trust region (if active)\r
991 #\r
992 switch $no_parm($opt,recompute_method) {\r
993 fd { $opt fdh }\r
994 fd2 { $opt fdh2 }\r
995 analytic { $opt h }\r
996 none {\r
997 if { $no_parm($opt,nsteps) > 1 } {\r
998 puts stdout "newopt: try prfo step-calculation again with previous hessian"\r
999 copy_object type=matrix from=$opt.ph to=$opt.h\r
1000 set restored_hessian 1\r
1001 }\r
1002 }\r
1003 ignore {\r
1004 puts stdout "newopt: take step regardless of prfo-failure"\r
1005 $opt cl 0.5*$no_parm($opt,maxlength)\r
1006 $opt s\r
1007 $opt fg\r
1008 switch $no_parm($opt,hessian_update) {\r
1009 bfgs { $opt bfu }\r
1010 powell { $opt pu }\r
1011 mursar { $opt msu }\r
1012 hybrid { $opt pmsu }\r
1013 fde { $opt efdh $no_parm($opt,var_list) }\r
1014 fde2 { $opt efdh2 $no_parm($opt,var_list) }\r
1015 default { chemerr "Baker: illegal hessian_update argument $no_parm($opt,hessian_update)" }\r
1016 }\r
1017 }\r
1018 fde {\r
1019 if { $no_parm($opt,var_list) != " " } {\r
1020 set eigenmode_list $no_parm($opt,var_list)\r
1021 } else {\r
1022 set eigenmode_list ""\r
1023 for { set i 1 } { $i <= [expr $no_inh($opt,neg) + $no_inh($opt,zero)] } { incr i 1 } {\r
1024 lappend eigenmode_list $i\r
1025 }\r
1026 }\r
1027 puts stdout "newopt: finite differencing of selected eigenmodes to improve the hessian"\r
1028 $opt efdh $eigenmode_list\r
1029 }\r
1030 fde2 {\r
1031 if { $no_parm($opt,var_list) != " " } {\r
1032 set eigenmode_list $no_parm($opt,var_list)\r
1033 } else {\r
1034 set eigenmode_list ""\r
1035 for { set i 1 } { $i <= [expr $no_inh($opt,neg) + $no_inh($opt,zero)] } { incr i 1 } {\r
1036 lappend eigenmode_list $i\r
1037 }\r
1038 }\r
1039 puts stdout "newopt: finite differencing of selected eigenmodes to improve the hessian"\r
1040 $opt efdh2 $eigenmode_list\r
1041 }\r
1042 default { chemerr " Baker: illegal recompute_method argument" }\r
1043 }\r
1044 }\r
1045 # diagonalise hessian\r
1046 $opt dh\r
1047 # check hessian character\r
1048 $opt hc\r
1049 #\r
1050 switch $no_parm($opt,trust_region) {\r
1051 active {\r
1052 set no_parm($opt,maxlength) [ expr 0.5 * $no_parm($opt,maxlength) ]\r
1053 if { $no_parm($opt,maxlength) < $no_parm($opt,stepmin) } {\r
1054 set no_parm($opt,maxlength) $no_parm($opt,stepmin)\r
1055 }\r
1056 }\r
1057 default { }\r
1058 }\r
1059 }\r
1060 2 {\r
1061 # the p-rfo step is definitely wrong, abort the job\r
1062 $opt log\r
1063 return -2\r
1064 }\r
1065 default {\r
1066 puts stdout $a\r
1067 chemerr "Baker: error in return code for prfo -${a}-"\r
1068 }\r
1069 }\r
1070 #\r
1071 # set hessian inertia and energy values\r
1072 #\r
1073 set no_inh($opt,pos_before) $no_inh($opt,pos)\r
1074 set no_parm($opt,old_energy) $no_parm($opt,new_energy)\r
1075 #\r
1076 newopt_save_restart_files $opt x f g h diis\r
1077 }\r
1078 #\r
1079 # save final hessian if required\r
1080 #\r
1081 switch $no_parm($opt,final_hessian) {\r
1082 undefined { }\r
1083 default { copy_matrix from=$opt.h to=$no_parm($opt,final_hessian) }\r
1084 }\r
1085 \r
1086 newopt_save_restart_files $opt x g f h diis\r
1087 return 0\r
1088 }\r
1089 \r
