1 |
#define EXTERN extern |
2 |
|
3 |
#include "pcwin.h" |
4 |
#include "pcmod.h" |
5 |
|
6 |
#include "energies.h" |
7 |
#include "derivs.h" |
8 |
#include "utility.h" |
9 |
#include "pot.h" |
10 |
#include "gmmx.h" |
11 |
|
12 |
// mode |
13 |
#define NONE 0 |
14 |
#define NEWTON 1 |
15 |
#define TNCG 2 |
16 |
#define DTNCG 3 |
17 |
#define AUTO 4 |
18 |
// method |
19 |
#define DIAG 5 |
20 |
#define SSOR 6 |
21 |
#define ICCG 7 |
22 |
#define BLOCK 8 |
23 |
|
24 |
|
25 |
#define Success 0 |
26 |
#define ReSearch 1 |
27 |
#define WideAngle 2 |
28 |
#define BadIntpln 3 |
29 |
#define IntplnErr 4 |
30 |
#define blank 5 |
31 |
#define Failure 6 |
32 |
|
33 |
void bounds(void); |
34 |
double energy(void); |
35 |
void tncg(int,int,int *,double *,double *, double, double (*)(), void (*)()); |
36 |
void mqn(int , int, int *,double *,double *, double *, double (*)() ); |
37 |
double minimize1(double *, double *); |
38 |
void search(int,double *,double *,double *,double *,double,double *,int *,double (*)(),int *); |
39 |
double newton1(double *, double *); |
40 |
void newton2(int, double *,double *,int *, int *, int *, double *); |
41 |
void hessian(int, double *, int *, int *, int *,double *); |
42 |
void gradient(void); |
43 |
void piseq(int, int); |
44 |
void piden(void); |
45 |
void inesc(char *); |
46 |
void outminstat(int , double ,double ); |
47 |
void RefreshScreen(void); |
48 |
void minimize(void); |
49 |
double minimiz1(double *, double *); |
50 |
|
51 |
struct t_minvar{ |
52 |
double cappa, stpmin, stpmax, angmax; |
53 |
int intmax; |
54 |
} minvar; |
55 |
|
56 |
EXTERN struct t_minim_control { |
57 |
int type, method, field, added_const; |
58 |
char added_path[256],added_name[256]; |
59 |
} minim_control; |
60 |
|
61 |
EXTERN struct t_minim_values { |
62 |
int iprint, ndc, nconst; |
63 |
float dielc; |
64 |
} minim_values; |
65 |
|
66 |
double scale2; |
67 |
|
68 |
void minimize() |
69 |
{ |
70 |
int i,nvar, iter, icount; |
71 |
double minimum,grdmin; |
72 |
double minimiz1(), newton1(); |
73 |
double etot; |
74 |
double *xx; // xx[maxvar]; |
75 |
int method, maxvar; |
76 |
|
77 |
maxvar = 3*natom; |
78 |
xx = dvector(0,maxvar); |
79 |
|
80 |
grdmin = 1.0; |
81 |
scale2 = 12.0; // bfgs 12.0 |
82 |
|
83 |
if (minim_control.method == 1 || minim_control.method == 3 || minim_control.method == 4) |
84 |
{ |
85 |
nvar = 0; |
86 |
icount = 0; |
87 |
for (i=1; i <= natom; i++) |
88 |
{ |
89 |
if (atom[i].use) |
90 |
{ |
91 |
xx[nvar] = atom[i].x*scale2; |
92 |
nvar++; |
93 |
xx[nvar] = atom[i].y*scale2; |
94 |
nvar++; |
95 |
xx[nvar] = atom[i].z*scale2; |
96 |
nvar++; |
97 |
} |
98 |
} |
99 |
|
100 |
method = 1; |
101 |
grdmin = 0.5; |
102 |
if (minim_control.method == 1) |
103 |
grdmin = 0.1; |
104 |
else |
105 |
grdmin = 0.5; |
106 |
|
107 |
mqn(nvar,method, &iter,xx, &minimum, &grdmin, minimiz1 ); |
108 |
|
109 |
} |
110 |
if (grdmin > 1.00) |
111 |
{ |
112 |
energies.total = 10000.0; |
113 |
free_dvector(xx, 0, maxvar); |
114 |
return; |
115 |
} |
116 |
if (minim_control.method == 2 || minim_control.method == 3 || minim_control.method == 4) |
117 |
{ |
118 |
scale2 = 1.0; // tcng |
119 |
grdmin = 0.0001; |
120 |
nvar = 0; |
121 |
icount = 0; |
122 |
for (i=1; i <= natom; i++) |
123 |
{ |
124 |
if (atom[i].use) |
125 |
{ |
126 |
xx[nvar] = atom[i].x*scale2; |
127 |
nvar++; |
128 |
xx[nvar] = atom[i].y*scale2; |
129 |
nvar++; |
130 |
xx[nvar] = atom[i].z*scale2; |
131 |
nvar++; |
132 |
} |
133 |
} |
134 |
|
135 |
tncg(nvar,method,&iter, xx, &minimum, grdmin,newton1, newton2); |
136 |
|
137 |
|
138 |
nvar = 0; |
139 |
for (i=1; i <= natom; i++) |
140 |
{ |
141 |
if (atom[i].use) |
142 |
{ |
143 |
atom[i].x = xx[nvar]/scale2; |
144 |
nvar++; |
145 |
atom[i].y = xx[nvar]/scale2; |
146 |
nvar++; |
147 |
atom[i].z = xx[nvar]/scale2; |
148 |
nvar++; |
149 |
} |
150 |
} |
151 |
} |
152 |
|
153 |
if (minimum < -1000.0) |
154 |
{ |
155 |
|
156 |
etot = energy(); |
157 |
free_dvector(xx, 0, maxvar); |
158 |
return; |
159 |
} |
160 |
free_dvector(xx, 0, maxvar); |
161 |
} |
162 |
|
163 |
void mqn(int nvar,int method,int *iter, double *x, double *minimum, double *grdmin, double (*fgvalue) ()) |
164 |
{ |
165 |
int i,ncalls,nerror; |
166 |
int niter,period,nstart; |
167 |
double fast,slow,epsln,d1temp,d2temp; |
168 |
double f,f_old,f_new,f_move; |
169 |
double rms,beta,x_move,g_norm,g_rms; |
170 |
double gg,gg_old; |
171 |
double sg,dg,sd,dd,angle; |
172 |
double *g, *p; // g[maxvar]; |
173 |
double *x_old, *g_old; // x_old[maxvar],g_old[maxvar]; |
174 |
double *s, *d; // p[maxvar],s[maxvar],d[maxvar]; |
175 |
double fctmin; |
176 |
int restart, terminate; |
177 |
int maxiter, nextiter,status, maxvar; |
178 |
|
179 |
maxvar = 3*natom; |
180 |
x_old = dvector(0,maxvar); |
181 |
g_old = dvector(0,maxvar); |
182 |
g = dvector(0,maxvar); |
183 |
p = dvector(0,maxvar); |
184 |
s = dvector(0,maxvar); |
185 |
d = dvector(0,maxvar); |
186 |
|
187 |
ncalls = 0; |
188 |
rms = sqrt((float)nvar)/ sqrt(3.0); |
189 |
restart = TRUE; |
190 |
terminate = FALSE; |
191 |
status = blank; |
192 |
nerror = 0; |
193 |
|
194 |
fctmin = -10000.0; |
195 |
maxiter = 1000; |
196 |
nextiter = 1; |
197 |
fast = 0.5; |
198 |
slow = 0.0; |
199 |
epsln = 1.0e-16; |
200 |
if (nvar > 200) |
201 |
period = nvar; |
202 |
else |
203 |
period = 200; |
204 |
minvar.cappa = .1; |
205 |
minvar.stpmin = 1.0e-20; |
206 |
minvar.stpmax = 5.0; |
207 |
minvar.angmax = 100.0; |
208 |
minvar.intmax = 5; |
209 |
|
210 |
niter = nextiter -1; |
211 |
maxiter = niter + maxiter; |
212 |
ncalls = ncalls + 1; |
213 |
f = fgvalue(x, g); // get function and first deriv at original point |
214 |
g_norm = 0.0; |
215 |
for (i=0; i < nvar; i++) |
216 |
{ |
217 |
x_old[i] = x[i]; |
218 |
g_old[i] = g[i]; |
219 |
g_norm += g[i]*g[i]; |
220 |
} |
221 |
g_norm = sqrt(g_norm); |
222 |
f_move = 0.5*minvar.stpmax*g_norm; |
223 |
g_rms = g_norm*scale2/rms; |
224 |
|
225 |
|
226 |
if (niter > maxiter) |
227 |
terminate = TRUE; |
228 |
if (f < fctmin) |
229 |
terminate = TRUE; |
230 |
if (g_rms < *grdmin) |
231 |
terminate = TRUE; |
232 |
|
233 |
while ( ! terminate) |
234 |
{ |
235 |
niter++; |
236 |
status = blank; |
237 |
|
238 |
|
239 |
if (restart || method == 0) |
240 |
{ |
241 |
for (i=0; i < nvar; i++) |
242 |
p[i] = -g[i]; |
243 |
nstart = niter; |
244 |
restart = FALSE; |
245 |
} else if (method == 1) // BFGS method |
246 |
{ |
247 |
sg = 0.0; |
248 |
dg = 0.0; |
249 |
dd = 0.0; |
250 |
sd = 0.0; |
251 |
for (i=0; i < nvar; i++) |
252 |
{ |
253 |
sg += s[i]*g[i]; |
254 |
dg += d[i]*g[i]; |
255 |
dd += d[i]*d[i]; |
256 |
sd += s[i]*d[i]; |
257 |
} |
258 |
for (i=0; i < nvar; i++) |
259 |
{ |
260 |
d1temp = (d[i]*sg + s[i]*dg)/sd; |
261 |
d2temp = (1.0+dd/sd)*(s[i]*sg/sd); |
262 |
p[i] = -g[i] + d1temp - d2temp; |
263 |
} |
264 |
} else if (method == 2) // Fletcher Reeves |
265 |
{ |
266 |
gg = 0.0; |
267 |
gg_old = 0.0; |
268 |
for (i=0; i < nvar; i++) |
269 |
{ |
270 |
gg += g[i]*g[i]; |
271 |
gg_old += g_old[i]*g_old[i]; |
272 |
} |
273 |
beta = gg/gg_old; |
274 |
for (i=0; i < nvar; i++) |
275 |
p[i] = -g[i] + beta*p[i]; |
276 |
} else if (method == 3) // Polak Ribere |
277 |
{ |
278 |
dg = 0.0; |
279 |
gg_old = 0.0; |
280 |
for (i=0; i < nvar; i++) |
281 |
{ |
282 |
dg += d[i]*g[i]; |
283 |
gg_old += g_old[i]*g_old[i]; |
284 |
} |
285 |
beta = dg/gg_old; |
286 |
for (i=0; i < nvar; i++) |
287 |
p[i] = -g[i] + beta*p[i]; |
288 |
} |
289 |
|
290 |
// do a line search |
291 |
f_old = f; |
292 |
search(nvar,&f,g,x,p,f_move,&angle,&ncalls,fgvalue,&status); |
293 |
if (status == Failure) |
294 |
{ |
295 |
g_rms = 1000.0; |
296 |
terminate = TRUE; |
297 |
goto L_DONE; |
298 |
} |
299 |
f_new = f; |
300 |
|
301 |
f_move = f_old - f_new; |
302 |
x_move = 0.0; |
303 |
g_norm = 0.0; |
304 |
for (i=0; i < nvar; i++) |
305 |
{ |
306 |
s[i] = x[i] - x_old[i]; |
307 |
d[i] = g[i] - g_old[i]; |
308 |
x_move += s[i]*s[i]; |
309 |
g_norm += g[i]*g[i]; |
310 |
x_old[i] = x[i]; |
311 |
g_old[i] = g[i]; |
312 |
} |
313 |
x_move = sqrt(x_move) / (scale2 * rms); |
314 |
g_norm = sqrt(g_norm); |
315 |
g_rms = g_norm * scale2/rms; |
316 |
|
317 |
// function increase |
318 |
if (f_move <= 0.0) |
319 |
{ |
320 |
// status = Increase; |
321 |
nerror = nerror + 1; |
322 |
if (nerror == 3) |
323 |
terminate = TRUE; |
324 |
else |
325 |
restart = TRUE; |
326 |
|
327 |
for(i=0; i < nvar; i++) |
328 |
{ |
329 |
x[i] = x_old[i]; |
330 |
g[i] = g_old[i]; |
331 |
} |
332 |
} |
333 |
if (x_move < epsln) |
334 |
{ |
335 |
nerror++; |
336 |
if (nerror > 3) |
337 |
terminate = TRUE; |
338 |
else |
339 |
restart = TRUE; |
340 |
} |
341 |
// normal termination |
342 |
if (f < fctmin) |
343 |
{ |
344 |
// status = SmallFct; |
345 |
terminate = TRUE; |
346 |
} |
347 |
if (g_rms < *grdmin) |
348 |
{ |
349 |
// status = SmallGrad; |
350 |
nerror++; |
351 |
if (nerror > 1) |
352 |
terminate = TRUE; |
353 |
else |
354 |
restart = TRUE; |
355 |
} |
356 |
|
357 |
} |
358 |
L_DONE: |
359 |
*minimum = f; |
360 |
*grdmin = g_rms; |
361 |
*iter = niter; |
362 |
free_dvector(x_old ,0,maxvar); |
363 |
free_dvector(g_old,0,maxvar); |
364 |
free_dvector( g ,0,maxvar); |
365 |
free_dvector( p ,0,maxvar); |
366 |
free_dvector( s ,0,maxvar); |
367 |
free_dvector( d ,0,maxvar); |
368 |
} |
369 |
|
370 |
double minimiz1(double *xx, double *g) |
371 |
{ |
372 |
int i,nvar; |
373 |
double e_min; |
374 |
|
375 |
nvar = 0; |
376 |
for (i = 1; i <= natom; i++) |
377 |
{ |
378 |
if (atom[i].use) |
379 |
{ |
380 |
atom[i].x = xx[nvar]/scale2; |
381 |
nvar++; |
382 |
atom[i].y = xx[nvar]/scale2; |
383 |
nvar++; |
384 |
atom[i].z = xx[nvar]/scale2; |
385 |
nvar++; |
386 |
} |
387 |
} |
388 |
|
389 |
gradient(); |
390 |
e_min = energies.total; |
391 |
|
392 |
nvar = 0; |
393 |
for (i=1; i <= natom; i++) |
394 |
{ |
395 |
if (atom[i].use) |
396 |
{ |
397 |
xx[nvar] = atom[i].x*scale2; |
398 |
g[nvar] = deriv.d1[i][0]/scale2; |
399 |
nvar++; |
400 |
xx[nvar] = atom[i].y*scale2; |
401 |
g[nvar] = deriv.d1[i][1]/scale2; |
402 |
nvar++; |
403 |
xx[nvar] = atom[i].z*scale2; |
404 |
g[nvar] = deriv.d1[i][2]/scale2; |
405 |
nvar++; |
406 |
} |
407 |
} |
408 |
return(e_min); |
409 |
} |
410 |
|
411 |
double newton1(double *xx, double *g) |
412 |
{ |
413 |
int i, nvar; |
414 |
double e; |
415 |
|
416 |
nvar = 0; |
417 |
for (i=1; i <= natom; i++) |
418 |
{ |
419 |
if (atom[i].use) |
420 |
{ |
421 |
atom[i].x = xx[nvar]; |
422 |
nvar++; |
423 |
atom[i].y = xx[nvar]; |
424 |
nvar++; |
425 |
atom[i].z = xx[nvar]; |
426 |
nvar++; |
427 |
} |
428 |
} |
429 |
|
430 |
gradient(); |
431 |
e = energies.total; |
432 |
|
433 |
nvar = 0; |
434 |
for (i=1; i <= natom; i++) |
435 |
{ |
436 |
if (atom[i].use) |
437 |
{ |
438 |
xx[nvar] = atom[i].x; |
439 |
g[nvar] = deriv.d1[i][0]; |
440 |
nvar++; |
441 |
xx[nvar] = atom[i].y; |
442 |
g[nvar] = deriv.d1[i][1]; |
443 |
nvar++; |
444 |
xx[nvar] = atom[i].z; |
445 |
g[nvar] = deriv.d1[i][2]; |
446 |
nvar++; |
447 |
} |
448 |
} |
449 |
return(e); |
450 |
} |
451 |
|
452 |
void newton2(int mode, double *xx,double *h,int *hinit, |
453 |
int *hstop, int *hindex, double *hdiag) |
454 |
{ |
455 |
int i,j,k,nvar, maxvar, nuse, maxhess; |
456 |
int *hvar, *huse; // hvar[maxvar],huse[maxvar]; |
457 |
|
458 |
if (mode == NONE) |
459 |
return; |
460 |
|
461 |
maxvar = 3*natom; |
462 |
hvar = ivector(0,maxvar); |
463 |
huse = ivector(0,maxvar); |
464 |
|
465 |
nvar = 0; |
466 |
nuse = TRUE; |
467 |
for (i=1; i <= natom; i++) |
468 |
{ |
469 |
if (atom[i].use) |
470 |
{ |
471 |
atom[i].x = xx[nvar]; |
472 |
nvar++; |
473 |
atom[i].y = xx[nvar]; |
474 |
nvar++; |
475 |
atom[i].z = xx[nvar]; |
476 |
nvar++; |
477 |
} else |
478 |
nuse = FALSE; |
479 |
} |
480 |
|
481 |
if (natom < 300) |
482 |
maxhess = (3*natom*(3*natom-1))/2; |
483 |
else if (natom < 800) |
484 |
maxhess = (3*natom*(3*natom-1))/3; |
485 |
else |
486 |
maxhess = (3*natom*(3*natom-1))/20; |
487 |
|
488 |
hessian(maxhess, h,hinit,hstop,hindex,hdiag); |
489 |
|
490 |
nvar = 0; |
491 |
if (nuse == FALSE) |
492 |
{ |
493 |
for (i=1; i <= natom; i++) |
494 |
{ |
495 |
k = 3*(i-1); |
496 |
if (atom[i].use) |
497 |
{ |
498 |
for (j=0; j < 3; j++) |
499 |
{ |
500 |
hvar[nvar] = j+k; |
501 |
huse[j+k] = nvar; |
502 |
nvar++; |
503 |
} |
504 |
} else |
505 |
{ |
506 |
for (j=0; j < 3; j++) |
507 |
huse[j+k] = 0; |
508 |
} |
509 |
} |
510 |
for (i=0; i < nvar; i++) |
511 |
{ |
512 |
k = hvar[i]; |
513 |
hinit[i] = hinit[k]; |
514 |
hstop[i] = hstop[k]; |
515 |
hdiag[i] = hdiag[k]; |
516 |
for (j=hinit[i]; j < hstop[i]; j++) |
517 |
hindex[j] = huse[hindex[j]]; |
518 |
} |
519 |
} |
520 |
// |
521 |
nvar = 0; |
522 |
for (i=1; i <= natom; i++) |
523 |
{ |
524 |
if (atom[i].use) |
525 |
{ |
526 |
xx[nvar] = atom[i].x; |
527 |
nvar++; |
528 |
xx[nvar] = atom[i].y; |
529 |
nvar++; |
530 |
xx[nvar] = atom[i].z; |
531 |
nvar++; |
532 |
} |
533 |
} |
534 |
free_ivector(hvar ,0,maxvar); |
535 |
free_ivector(huse ,0,maxvar); |
536 |
|
537 |
} |
538 |
|
539 |
|
540 |
|