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, 4 months ago) by gilbertke
File size: 12062 byte(s)
Log Message:
major rewrite - removing global data, adding electrostatics tag to read_sdf
Line File contents
1 #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