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