ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/branches/sync4pymol12/src/mengine/src/minimize.c
Revision: 121
Committed: Tue Jun 16 22:35:06 2009 UTC (12 years, 11 months ago) by wdelano
File size: 13023 byte(s)
Log Message:
updated mengine
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     struct timeval tv;
57 wdelano 112
58     maxvar = 3*natom;
59     xx = dvector(0,maxvar);
60    
61     grdmin = 1.0;
62     scale2 = 12.0; // bfgs 12.0
63    
64 wdelano 121 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 wdelano 112 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 wdelano 121 mqn(start_time,nvar,method, &iter,xx, &minimum, &grdmin, minimiz1 );
96 wdelano 112 }
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 wdelano 121 tncg(start_time, nvar,method,&iter, xx, &minimum, grdmin,minimiz1, newton2);
122 wdelano 112
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 wdelano 121 if (job_control.monitor)
148     job_control.monitor = FALSE;
149 wdelano 112 }
150     // ======================================
151 wdelano 121 void mqn(double start_time,int nvar,int method,int *iter, double *x, double *minimum, double *grdmin, double (*fgvalue) ())
152 wdelano 112 {
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 wdelano 121 struct timeval tv;
167     double end_time,old_time,ttime;
168 wdelano 112
169 wdelano 121 old_time = start_time;
170 wdelano 112 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 wdelano 121 // 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 wdelano 112
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