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