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