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, 11 months ago) by wdelano
File size: 13358 byte(s)
Log Message:
Line User Rev File contents
1 wdelano 112 #define EXTERN extern
2    
3     #include "pcwin.h"
4     #include "pcmod.h"
5     #include "utility.h"
6 wdelano 121 #include "job_control.h"
7     #include <sys/timeb.h>
8     #include <time.h>
9 wdelano 112
10     // mode
11     #define NONE 0
12     #define Failure 6
13    
14     double energy(void);
15 wdelano 121 void tncg(double,int,int,int *,double *,double *, double, double (*)(), void (*)());
16     void mqn(double,int , int, int *,double *,double *, double *, double (*)() );
17 wdelano 112 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 wdelano 121 void write_sdf(int);
28 wdelano 112
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 wdelano 121 // double minimiz1();
52     double etot,stime;
53 wdelano 112 double *xx; // xx[maxvar];
54     int method, maxvar;
55 wdelano 121 double start_time;
56 wdelano 112
57     maxvar = 3*natom;
58     xx = dvector(0,maxvar);
59    
60     grdmin = 1.0;
61     scale2 = 12.0; // bfgs 12.0
62    
63 wdelano 121 if (job_control.monitor)
64     {
65 wdelano 122 #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 wdelano 121 }
75     else
76     start_time = 0;
77 wdelano 112 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 wdelano 121 mqn(start_time,nvar,method, &iter,xx, &minimum, &grdmin, minimiz1 );
102 wdelano 112 }
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 wdelano 121 tncg(start_time, nvar,method,&iter, xx, &minimum, grdmin,minimiz1, newton2);
128 wdelano 112
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 wdelano 121 if (job_control.monitor)
154     job_control.monitor = FALSE;
155 wdelano 112 }
156     // ======================================
157 wdelano 121 void mqn(double start_time,int nvar,int method,int *iter, double *x, double *minimum, double *grdmin, double (*fgvalue) ())
158 wdelano 112 {
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 wdelano 121 double end_time,old_time,ttime;
173 wdelano 112
174 wdelano 121 old_time = start_time;
175 wdelano 112 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 wdelano 121 // check monitor
353     if (job_control.monitor)
354     {
355 wdelano 122 #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 wdelano 121 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 wdelano 112
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