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