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 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 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 status = 0;
168 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 status = 0;
213
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