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