mengine/src/minimize.c

#define EXTERN extern

#include "pcwin.h"
#include "pcmod.h"

#include "energies.h"
#include "derivs.h"
#include "utility.h"
#include "pot.h"
#include "gmmx.h"

// mode
#define NONE      0
#define  Failure     6

double energy(void);
void tncg(int,int,int *,double *,double *, double, double (*)(),  void (*)());
void mqn(int , int, int *,double *,double *, double *, double (*)() );
double minimize1(double *, double *);
void search(int,double *,double *,double *,double *,double,double *,int *,double (*)(),int *);
double newton1(double *, double *);
void newton2(int, double *,double *,int *, int *, int *, double *);
void hessian(int, double *, int *, int *, int *,double *);
void gradient(void);
void minimize(void);
double minimiz1(double *, double *);

struct t_minvar{
     double cappa, stpmin, stpmax, angmax;
     int   intmax;
     }  minvar;
     
EXTERN struct t_minim_control {
        int type, method, field, added_const;
        char added_path[256],added_name[256];
        } minim_control;
        
EXTERN struct t_minim_values {
        int iprint, ndc, nconst;
        float dielc;
        } minim_values;
    
double scale2;

void minimize()
{
   int i,nvar, iter, icount;
   double minimum,grdmin;
   double minimiz1(), newton1();
   double etot;
   double *xx;  //  xx[maxvar];
   int method, maxvar;
   
   maxvar = 3*natom;
   xx = dvector(0,maxvar);
      
   grdmin = 1.0;
   scale2 = 12.0; // bfgs  12.0

   if (minim_control.method == 1 || minim_control.method == 3 || minim_control.method == 4)
   {
     nvar = 0;
     icount = 0;
     for (i=1; i <= natom; i++)
     {
       if (atom[i].use)
       {
         xx[nvar] = atom[i].x*scale2;
         nvar++;
         xx[nvar] = atom[i].y*scale2;
         nvar++;
         xx[nvar] = atom[i].z*scale2;
         nvar++;
       }
     }

     method = 1;
     grdmin = 0.5;
     if (minim_control.method == 1)
        grdmin = 0.1;
     else
        grdmin = 0.5;
       
     mqn(nvar,method, &iter,xx, &minimum, &grdmin, minimiz1 );
    
   }
   if (grdmin > 1.00)
   {
       energies.total = 10000.0;
       free_dvector(xx, 0, maxvar);
       return;
   }
      if (minim_control.method == 2 || minim_control.method == 3 || minim_control.method == 4)
      {
        scale2 = 1.0;   // tcng
        grdmin = 0.0001;
        nvar = 0;
        icount = 0;
        for (i=1; i <= natom; i++)
        {
          if (atom[i].use)
          {
            xx[nvar] = atom[i].x*scale2;
            nvar++;
            xx[nvar] = atom[i].y*scale2;
            nvar++;
            xx[nvar] = atom[i].z*scale2;
            nvar++;
          }
        }
        
        tncg(nvar,method,&iter, xx, &minimum, grdmin,newton1, newton2); 
       
   
        nvar = 0;
        for (i=1; i <= natom; i++)
        {
          if (atom[i].use)
          {
            atom[i].x =  xx[nvar]/scale2;
            nvar++;
            atom[i].y =  xx[nvar]/scale2;
            nvar++;
            atom[i].z =  xx[nvar]/scale2;
            nvar++;
          }
        }
      }
   
   if (minimum < -1000.0)
   {
    
       etot = energy();    
       free_dvector(xx, 0, maxvar);
       return;      
   }   
   free_dvector(xx, 0, maxvar);
}

void mqn(int nvar,int method,int *iter, double *x,  double *minimum, double *grdmin, double (*fgvalue) ())
{
      int i,ncalls,nerror;
      int niter,period,nstart;
      double fast,slow,epsln,d1temp,d2temp;
      double f,f_old,f_new,f_move;
      double rms,beta,x_move,g_norm,g_rms;
      double gg,gg_old;
      double sg,dg,sd,dd,angle;
      double *g, *p;                         //   g[maxvar];
      double *x_old, *g_old;             //   x_old[maxvar],g_old[maxvar];
      double  *s, *d;                 //   p[maxvar],s[maxvar],d[maxvar];
      double fctmin;
      int restart, terminate;
      int maxiter, nextiter,status, maxvar;

      maxvar = 3*natom;
      x_old = dvector(0,maxvar);
      g_old = dvector(0,maxvar);
      g = dvector(0,maxvar);
      p = dvector(0,maxvar);
      s = dvector(0,maxvar);
      d = dvector(0,maxvar);
      
      ncalls = 0;
      rms = sqrt((float)nvar)/ sqrt(3.0);
      restart = TRUE;
      terminate = FALSE;
      status = 0;
      nerror = 0;

      fctmin = -10000.0;
      maxiter = 1000;
      nextiter = 1;
      fast = 0.5;
      slow = 0.0;
      epsln = 1.0e-16;
      if (nvar > 200)
         period = nvar;
      else
         period = 200;
      minvar.cappa = .1;
      minvar.stpmin = 1.0e-20;
      minvar.stpmax = 5.0;
      minvar.angmax = 100.0;
      minvar.intmax = 5;

      niter = nextiter -1;
      maxiter = niter + maxiter;
      ncalls = ncalls + 1;
      f = fgvalue(x, g); // get function and first deriv at original point
      g_norm = 0.0;
      for (i=0; i < nvar; i++)
      {
          x_old[i] = x[i];
          g_old[i] = g[i];
          g_norm += g[i]*g[i];
      }
      g_norm = sqrt(g_norm);
      f_move = 0.5*minvar.stpmax*g_norm;
      g_rms = g_norm*scale2/rms;


     if (niter > maxiter)
        terminate = TRUE;
     if (f < fctmin)
        terminate = TRUE;
     if (g_rms < *grdmin)
         terminate = TRUE;

     while ( ! terminate)
     {
         niter++;
         status = 0;   


         if (restart || method == 0)
         {
             for (i=0; i < nvar; i++)
                 p[i] = -g[i];
             nstart = niter;
             restart = FALSE;
         } else if (method == 1) // BFGS method
         {
             sg = 0.0;
             dg = 0.0;
             dd = 0.0;
             sd = 0.0;
             for (i=0; i < nvar; i++)
             {
               sg += s[i]*g[i];
               dg += d[i]*g[i];
               dd += d[i]*d[i];
               sd += s[i]*d[i];
             }
             for (i=0; i < nvar; i++)
             {
                 d1temp = (d[i]*sg + s[i]*dg)/sd;
                 d2temp = (1.0+dd/sd)*(s[i]*sg/sd);
                 p[i] = -g[i] + d1temp - d2temp;
             }
         } else if (method == 2)  // Fletcher Reeves
         {
             gg = 0.0;
             gg_old = 0.0;
             for (i=0; i < nvar; i++)
             {
                 gg += g[i]*g[i];
                 gg_old += g_old[i]*g_old[i];
             }
             beta = gg/gg_old;
             for (i=0; i < nvar; i++)
                p[i] = -g[i] + beta*p[i];
         } else if (method == 3)  // Polak Ribere
         {
             dg = 0.0;
             gg_old = 0.0;
             for (i=0; i < nvar; i++)
             {
                 dg += d[i]*g[i];
                 gg_old += g_old[i]*g_old[i];
             }
             beta = dg/gg_old;
             for (i=0; i < nvar; i++)
                p[i] = -g[i] + beta*p[i];
        }

         // do a line search
         f_old = f;
         search(nvar,&f,g,x,p,f_move,&angle,&ncalls,fgvalue,&status);
         if (status == Failure)
         {
             g_rms = 1000.0;
             terminate = TRUE;
             goto L_DONE;
         }
         f_new = f;
         
         f_move = f_old - f_new;
         x_move = 0.0;
         g_norm = 0.0;
         for (i=0; i < nvar; i++)
         {
            s[i] = x[i] - x_old[i];
            d[i] = g[i] - g_old[i];
            x_move += s[i]*s[i];
            g_norm += g[i]*g[i];
            x_old[i] = x[i];
            g_old[i] = g[i];
         }
         x_move = sqrt(x_move) / (scale2 * rms);
         g_norm = sqrt(g_norm);
         g_rms = g_norm * scale2/rms;
         
// function increase
         if (f_move <= 0.0)
         {
 //           status = Increase;
            nerror = nerror + 1;
            if (nerror == 3)
               terminate = TRUE;
            else
               restart = TRUE;

            for(i=0; i < nvar; i++)
            {
               x[i] = x_old[i];
               g[i] = g_old[i];
            }
         }
         if (x_move < epsln)
         {
             nerror++;
             if (nerror > 3)
               terminate = TRUE;
             else
               restart = TRUE;
         }
// normal termination
         if (f < fctmin)
         {
 //           status = SmallFct;
            terminate = TRUE;
         }
         if (g_rms < *grdmin)
         {
//            status = SmallGrad;
            nerror++;
            if (nerror > 1)
             terminate = TRUE;
            else
             restart = TRUE;
         }

     }
L_DONE:
     *minimum = f;
     *grdmin = g_rms;
     *iter = niter;
     free_dvector(x_old ,0,maxvar);
     free_dvector(g_old,0,maxvar);
     free_dvector( g ,0,maxvar);
     free_dvector( p ,0,maxvar);
     free_dvector( s ,0,maxvar);
     free_dvector( d ,0,maxvar);             
}

double minimiz1(double *xx, double *g)
{
    int i,nvar;
    double e_min;

    nvar = 0;
    for (i = 1; i <= natom; i++)
    {
        if (atom[i].use)
        {
            atom[i].x = xx[nvar]/scale2;
            nvar++;
            atom[i].y = xx[nvar]/scale2;
            nvar++;
            atom[i].z = xx[nvar]/scale2;
            nvar++;
        }
    }

    gradient();
    e_min = energies.total;

    nvar = 0;
    for (i=1; i <= natom; i++)
    {
        if (atom[i].use)
        {
            xx[nvar] = atom[i].x*scale2;
            g[nvar] = deriv.d1[i][0]/scale2;
            nvar++;
            xx[nvar] = atom[i].y*scale2;
            g[nvar] = deriv.d1[i][1]/scale2;
            nvar++;
            xx[nvar] = atom[i].z*scale2;
            g[nvar] = deriv.d1[i][2]/scale2;
            nvar++;
        }
    }
    return(e_min);
}

double newton1(double *xx, double *g)
{
    int i, nvar;
    double e;

    nvar = 0;
    for (i=1; i <= natom; i++)
    {
        if (atom[i].use)
        {
            atom[i].x = xx[nvar];
            nvar++;
            atom[i].y = xx[nvar];
            nvar++;
            atom[i].z = xx[nvar];
            nvar++;
        }
    }

    gradient();
    e = energies.total;

    nvar = 0;
    for (i=1; i <= natom; i++)
    {
        if (atom[i].use)
        {
            xx[nvar] = atom[i].x;
            g[nvar] = deriv.d1[i][0];
            nvar++;
            xx[nvar] = atom[i].y;
            g[nvar] = deriv.d1[i][1];
            nvar++;
            xx[nvar] = atom[i].z;
            g[nvar] = deriv.d1[i][2];
            nvar++;
        }
    }
    return(e);
}

void newton2(int mode, double *xx,double *h,int *hinit,
               int *hstop, int *hindex, double *hdiag)
{
    int i,j,k,nvar, maxvar, nuse, maxhess;
    int *hvar, *huse;   //  hvar[maxvar],huse[maxvar];

    if (mode == NONE)
       return;

    maxvar = 3*natom;
    hvar = ivector(0,maxvar);
    huse = ivector(0,maxvar);
    
    nvar = 0;
    nuse = TRUE;
    for (i=1; i <= natom; i++)
    {
        if (atom[i].use)
        {
           atom[i].x = xx[nvar];
           nvar++;
           atom[i].y = xx[nvar];
           nvar++;
           atom[i].z = xx[nvar];
           nvar++;
        } else
           nuse = FALSE;
    }

    if (natom < 300)
      maxhess = (3*natom*(3*natom-1))/2;
    else if (natom < 800)
      maxhess = (3*natom*(3*natom-1))/3;
    else
      maxhess = (3*natom*(3*natom-1))/20;

    hessian(maxhess, h,hinit,hstop,hindex,hdiag);

    nvar = 0;
    if (nuse == FALSE)
    {
       for (i=1; i <= natom; i++)
       {
           k = 3*(i-1);
           if (atom[i].use)
           {
              for (j=0; j < 3; j++)
              {
                hvar[nvar] = j+k;
                huse[j+k] = nvar;
                nvar++;
              }
           } else
           {
              for (j=0; j < 3; j++)
                 huse[j+k] = 0;
           }
       }
       for (i=0; i < nvar; i++)
       {
          k = hvar[i];
          hinit[i] = hinit[k];
          hstop[i] = hstop[k];
          hdiag[i] = hdiag[k];
          for (j=hinit[i]; j < hstop[i]; j++)
            hindex[j] = huse[hindex[j]];
       }
    }
//
    nvar = 0;
    for (i=1; i <= natom; i++)
    {
        if (atom[i].use)
        {
            xx[nvar] = atom[i].x;
            nvar++;
            xx[nvar] = atom[i].y;
            nvar++;
            xx[nvar] = atom[i].z;
            nvar++;
        }
    }         
    free_ivector(hvar ,0,maxvar);
    free_ivector(huse ,0,maxvar);

}
            
            
Revision:	63
Committed:	Wed Dec 3 03:46:32 2008 UTC (13 years, 7 months ago) by gilbertke
File size:	12212 byte(s)
Log Message:	updated read_sdf type_mmx and first pass at using best practices
Line	User	Rev	File contents
1	tjod	3	#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	gilbertke	63	status = 0;
169	tjod	3	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.5minvar.stpmaxg_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	gilbertke	63	status = 0;
214	tjod	3
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 = (3natom(3*natom-1))/2;
460			else if (natom < 800)
461			maxhess = (3natom(3*natom-1))/3;
462			else
463			maxhess = (3natom(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