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