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, 10 months ago) by gilbertke
File size: 12212 byte(s)
Log Message:
updated read_sdf type_mmx and first pass at using best practices
Line File contents
1 #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 status = 0;
169 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 status = 0;
214
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