ViewVC Help
View File | Revision Log | Show Annotations | View Changeset | Root Listing
root/freemol/trunk/src/mengine/src/minimize.c
Revision: 27
Committed: Tue Jul 8 19:10:25 2008 UTC (12 years, 4 months ago) by tjod
File size: 12675 byte(s)
Log Message:
move mengine src to reflect new freemol directory structure

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