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

Diff Legend

Removed lines
+ Added lines
< Changed lines
> Changed lines