ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/smi23d/src/mengine/minimize.c
Revision: 3
Committed: Mon Jun 9 21:38:26 2008 UTC (11 years, 4 months ago) by tjod
File size: 12675 byte(s)
Log Message:
test

Line File contents
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