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