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