//#include <mpi.h>
#include <cmath>
#include <iostream>
#include <fstream>
#include <iomanip>
#include "vectormatrixclass.h"

using namespace std;

//Step length for numerical derivation
const double h = .0001;
const double hh = h*h;
int antMC_global;
int Z_global;
double delta_t_global;
long idum_global;


static double sqrarg;
#define SQR(a) ((sqrarg=(a)) == 0.0 ? 0.0 : sqrarg*sqrarg)


static double maxarg1,maxarg2;
#define FMAX(a,b) (maxarg1=(a),maxarg2=(b),(maxarg1) > (maxarg2) ?\
        (maxarg1) : (maxarg2))


#define ALF 1.0e-4
#define TOLX 1.0e-7

void lnsrch(int n, Vector &xold, double fold, Vector &g, Vector &p, Vector &x,
	    double *f, double stpmax, int *check, double (*func)(Vector &p))
{
  int i;
  double a,alam,alam2,alamin,b,disc,f2,fold2,rhs1,rhs2,slope,sum,temp,
    test,tmplam;

  *check=0;
  for (sum=0.0,i=0;i<n;i++) sum += p(i)*p(i);
  sum=sqrt(sum);
  if (sum > stpmax)
    for (i=0;i<n;i++) p(i) *= stpmax/sum;
  for (slope=0.0,i=0;i<n;i++)
    slope += g(i)*p(i);
  test=0.0;
  for (i=0;i<n;i++) {
    temp=fabs(p(i))/FMAX(fabs(xold(i)),1.0);
    if (temp > test) test=temp;
  }
  alamin=TOLX/test;
  alam=1.0;
  for (;;) {
    for (i=0;i<n;i++) x(i)=xold(i)+alam*p(i);
    *f=(*func)(x);
    if (alam < alamin) {
      for (i=0;i<n;i++) x(i)=xold(i);
      *check=1;
      return;
    } else if (*f <= fold+ALF*alam*slope) return;
    else {
      if (alam == 1.0)
	tmplam = -slope/(2.0*(*f-fold-slope));
      else {
	rhs1 = *f-fold-alam*slope;
	rhs2=f2-fold2-alam2*slope;
	a=(rhs1/(alam*alam)-rhs2/(alam2*alam2))/(alam-alam2);
	b=(-alam2*rhs1/(alam*alam)+alam*rhs2/(alam2*alam2))/(alam-alam2);
	if (a == 0.0) tmplam = -slope/(2.0*b);
	else {
	  disc=b*b-3.0*a*slope;
	  if (disc<0.0) cout << "Roundoff problem in lnsrch." << endl;
	  else tmplam=(-b+sqrt(disc))/(3.0*a);
	}
	if (tmplam>0.5*alam)
	  tmplam=0.5*alam;
      }
    }
    alam2=alam;
    f2 = *f;
    fold2=fold;
    alam=FMAX(tmplam,0.1*alam);
  }
}
#undef ALF
#undef TOLX


#define ITMAX 200
#define EPS 3.0e-8
#define TOLX (4*EPS)
#define STPMX 100.0


void dfpmin(Vector &p, int n, double gtol, int *iter, double *fret,
	    double(*func)(Vector &p), void (*dfunc)(Vector &p, Vector &g))
{

  int check,i,its,j;
  double den,fac,fad,fae,fp,stpmax,sum=0.0,sumdg,sumxi,temp,test;
  Vector dg(n), g(n), hdg(n), pnew(n), xi(n);
  Matrix hessian(n,n);

  fp=(*func)(p);
  (*dfunc)(p,g);
  for (i = 0;i < n;i++) {
    for (j = 0; j< n;j++) hessian(i,j)=0.0;
    hessian(i,i)=1.0;
    xi(i) = -g(i);
    sum += p(i)*p(i);
  }
  stpmax=STPMX*FMAX(sqrt(sum),(double)n);
  for (its=1;its<=ITMAX;its++) {
    *iter=its;
    lnsrch(n,p,fp,g,xi,pnew,fret,stpmax,&check,func);
    fp = *fret;
    for (i = 0; i< n;i++) {
      xi(i)=pnew(i)-p(i);
      p(i)=pnew(i);
    }
    test=0.0;
    for (i = 0;i< n;i++) {
      temp=fabs(xi(i))/FMAX(fabs(p(i)),1.0);
      if (temp > test) test=temp;
    }
    if (test < TOLX) {
      return;
    }
    for (i=0;i<n;i++) dg(i)=g(i);
    (*dfunc)(p,g);
    test=0.0;
    den=FMAX(*fret,1.0);
    for (i=0;i<n;i++) {
      temp=fabs(g(i))*FMAX(fabs(p(i)),1.0)/den;
      if (temp > test) test=temp;
    }
    if (test < gtol) {
      return;
    }
    for (i=0;i<n;i++) dg(i)=g(i)-dg(i);
    for (i=0;i<n;i++) {
      hdg(i)=0.0;
      for (j=0;j<n;j++) hdg(i) += hessian(i,j)*dg(j);
    }
    fac=fae=sumdg=sumxi=0.0;
    for (i=0;i<n;i++) {
      fac += dg(i)*xi(i);
      fae += dg(i)*hdg(i);
      sumdg += SQR(dg(i));
      sumxi += SQR(xi(i));
    }
    if (fac*fac > EPS*sumdg*sumxi) {
      fac=1.0/fac;
      fad=1.0/fae;
      for (i=0;i<n;i++) dg(i)=fac*xi(i)-fad*hdg(i);
      for (i=0;i<n;i++) {
	for (j=0;j<n;j++) {
	  hessian(i,j) += fac*xi(i)*xi(j)
	    -fad*hdg(i)*hdg(j)+fae*dg(i)*dg(j);
	}
      }
    }
    for (i=0;i<n;i++) {
      xi(i)=0.0;
      for (j=0;j<n;j++) xi(i) -= hessian(i,j)*g(j);
    }
  }
  cout << "too many iterations in dfpmin" << endl;
}
#undef ITMAX
#undef EPS
#undef TOLX
#undef STPMX

     /*
     ** The function
     **           ran1()
     ** is an "Minimal" random number generator of Park and Miller
     ** (see Numerical recipe page 280) with Bays-Durham shuffle and
     ** added safeguards. Call with idum a negative integer to initialize;
     ** thereafter, do not alter idum between sucessive deviates in a
     ** sequence. RNMX should approximate the largest floating point value
     ** that is less than 1.
     ** The function returns a uniform deviate between 0.0 and 1.0
     ** (exclusive of end-point values).
     */

#define IA 16807
#define IM 2147483647
#define AM (1.0/IM)
#define IQ 127773
#define IR 2836
#define NTAB 32
#define NDIV (1+(IM-1)/NTAB)
#define EPS 1.2e-7
#define RNMX (1.0-EPS)

double ran1(long *idum)
{
   int             j;
   long            k;
   static long     iy=0;
   static long     iv[NTAB];
   double          temp;

   if (*idum <= 0 || !iy) {
      if (-(*idum) < 1) *idum=1;
      else *idum = -(*idum);
      for(j = NTAB + 7; j >= 0; j--) {
         k     = (*idum)/IQ;
         *idum = IA*(*idum - k*IQ) - IR*k;
         if(*idum < 0) *idum += IM;
         if(j < NTAB) iv[j] = *idum;
      }
      iy = iv[0];
   }
   k     = (*idum)/IQ;
   *idum = IA*(*idum - k*IQ) - IR*k;
   if(*idum < 0) *idum += IM;
   j     = iy/NDIV;
   iy    = iv[j];
   iv[j] = *idum;
   if((temp=AM*iy) > RNMX) return RNMX;
   else return temp;
}
#undef IA
#undef IM
#undef AM
#undef IQ
#undef IR
#undef NTAB
#undef NDIV
#undef EPS
#undef RNMX

// End: function ran1()





//Gaussian random deviate
double gaussian_deviate(long* idum)
{
  static int iset = 0;
  static double gset;
  double fac, rsq, v1, v2;

  if ( idum < 0) iset =0;
  if (iset == 0) {
    do {
      v1 = 2.*ran1(idum) -1.0;
      v2 = 2.*ran1(idum) -1.0;
      rsq = v1*v1+v2*v2;
    } while (rsq >= 1.0 || rsq == 0.);
    fac = sqrt(-2.*log(rsq)/rsq);
    gset = v1*fac;
    iset = 1;
    return v2*fac;
  } else {
    iset =0;
    return gset;
  }
}

//Distance from nucleus and electron i
double r_i(Matrix &R, int i) {
  
  double result = 0;
  for (int j = 0; j < 3; j++)
    result += R(i,j)*R(i,j);
  
  return sqrt(result);
}

//Distance between r_i and r_j
double r_ij(Matrix &R, int i, int j) {

  double result = 0;
  for (int k = 0; k < 3; k++)
    result += (R(i,k)-R(j,k)) * ((R(i,k)-R(j,k)));

  return sqrt(result);
}


//Hydrogen wave functions
//phi_j(r_i)  
double phi(int j, Matrix &R, int i, double alpha) {

  if (j == 0) {
    return exp(-alpha*r_i(R,i));
  }

  if (j == 1) {
    double r = r_i(R,i);
    return (1-alpha*r/2)*exp(-alpha*r/2);
  }

  if((j > 1) && (j < 5)) {
    return alpha*R(i,j-2)*exp(-alpha*r_i(R,i)/2);
  }

  cout << "ERROR: undefined j in phi" << endl;
}

//Derivative Hydrogen wave function
//dphi_j(r_i)/dr_ik
double phi_deriv(int j, Matrix &R, int i, int k, double alpha) {

  if (j == 0) {
    double r = r_i(R,i);
    return -alpha*R(i,k)*exp(-alpha*r) / r;
  }

  if (j == 1) {
    double r = r_i(R,i);
    return -alpha*R(i,k)*(2-alpha*r/2)*exp(-alpha*r/2) / (2*r);
  }

  if((j > 1) && (j < 5)) {
    if (k == j - 2) {
      double r = r_i(R,i);
      return alpha*exp(-alpha*r/2) * (1 - alpha*R(i,k)*R(i,k)/(2*r));
    }
    else {
      double r = r_i(R,i);
      return -alpha*alpha*R(i,j-2)*R(i,k)*exp(-alpha*r/2) / (2*r);
    }
  }

  cout << "ERROR: undefined j in phi_deriv" << endl;
}

//Second derivative of hydrogen-like wave function
//d^2phi_j(r_i)/dr_ik^2
double phi_deriv2(int j, Matrix &R, int i, int k, double alpha) {

  if (j == 0) {
    double r = r_i(R,i);
    return (pow(alpha*R(i,k)/r,2) - alpha*(r*r-R(i,k)*R(i,k))/(r*r*r))*exp(-alpha*r);
  }

  if (j == 1) {
    double r = r_i(R,i);
    return -alpha/2 * ( (r*r-R(i,k)*R(i,k))/(r*r*r)*(2-alpha*r/2)
			- alpha/2*pow(R(i,k)/r,2)*(3-alpha*r/2) )*exp(-alpha*r/2);
  }

  if((j > 1) && (j < 5)) {
    if (k == j - 2) {
      double r = r_i(R,i);
      return alpha*alpha*exp(-alpha*r/2)/r * ( pow(R(i,k),3)*(1/r+.5*alpha)/(2*r) - 1.5*R(i,k) );
    }
    else {
      double r = r_i(R,i);
      return alpha*alpha*R(i,j-2)*exp(-alpha*r/2)/2 * ( .5*alpha*pow(R(i,k)/r,2) - (r*r - R(i,k)*R(i,k))/pow(r,3) );
    }
  } 

  cout << "ERROR: undefined j in phi_deriv2" << endl;
}

//Determinant function
double determinant(Matrix &A, int dim) {

  if (dim == 2)
    return A(0,0)*A(1,1) - A(0,1)*A(1,0);

  double sum = 0;

  /*  for (int i = 0; i < dim; i++) {

    sub = new double*[dim-1];
    for (int j = 0; j < i; j++)
      sub[j] = &A[j][1];
    for (int j = i+1; j < dim; j++)
      sub[j-1] = &A[j][1];

    if(i % 2 == 0)
      sum += A[i][0] * determinant(sub,dim-1);
    else
      sum -= A[i][0] * determinant(sub,dim-1);

    delete[] sub;
  }
  return sum;
  */
}

//Slater-determinant
double slater(Matrix &R, double alpha, double Z, double Z2) {

  Matrix DUp(Z2,Z2), DDown(Z2,Z2);
  for (int i = 0; i < Z2; i++) {
    for (int j = 0; j < Z2; j++) {
      DUp(i,j) = phi(j,R,i,alpha);
      DDown(i,j) = phi(j,R,i+Z2,alpha);
    }
  }
  //Return product of spin up and spin down matrices
  double det = determinant(DUp,Z2)*determinant(DDown,Z2);
  return det;
}


void ludcmp(Matrix &a, int n, Vector &indx, double *d)
{
   int      i, imax, j, k;
   double   big, dum, sum, temp;
   Vector vv(n);

   *d = 1.0;                              // no row interchange yet
   for(i = 0; i < n; i++) {     // loop over rows to get scaling information
      big = 1.0E-10;
      for(j = 0; j < n; j++) {
	if((temp = fabs(a(i,j))) > big) big = temp;
      }
      if(big == 1.0E-10) {
         printf("\n\nSingular matrix in routine ludcmp()\n");
      }               
      vv(i) = 1.0/big;                 // save scaling */
   } // end i-loop */

   for(j = 0; j < n; j++) {     // loop over columns of Crout's method
      for(i = 0; i< j; i++) {   // not i = j
	sum = a(i,j);    
	for(k = 0; k < i; k++) sum -= a(i,k)*a(k,j);
	a(i,j) = sum;
      }
      big = 1.0E-10;   // initialization for search for largest pivot element
      for(i = j; i< n; i++) {
	sum = a(i,j);
	for(k = 0; k < j; k++) sum -= a(i,k)*a(k,j);
	a(i,j) = sum;
	if((dum = vv(i)*fabs(sum)) >= big) {
  	    big = dum;
	    imax = i;
	 }
      } // end i-loop
      if(j != imax) {    // do we need to interchange rows ?
         for(k = 0;k< n; k++) {       // yes
	   dum        = a(imax,k);
	   a(imax,k) = a(j,k);
	   a(j,k)    = dum;
	 }
	 (*d)    *= -1;            // and change the parit of d
	 vv(imax) = vv(j);         // also interchange scaling factor 
      }
      indx(j) = imax;
      if(fabs(a(j,j)) < 1.0E-10)  a(j,j) = 1.0E-10;

        /*
        ** if the pivot element is zero the matrix is singular
        ** (at least to the precision of the algorithm). For 
        ** some application of singular matrices, it is desirable
        ** to substitute ZERO for zero,
        */

      if(j < (n - 1)) {                   // divide by pivot element 
	dum = 1.0/a(j,j);
	for(i=j+1;i < n; i++) a(i,j) *= dum;
      }
   } // end j-loop over columns

}  // End: function ludcmp()

void lubksb(Matrix &a, int n, Vector &indx, Vector &b)
{
   int        i, ii = -1, ip, j;
   double     sum;

   for(i = 0; i< n; i++) {
     ip    = indx(i);
     sum   = b(ip);
     b(ip) = b(i);
     if(ii > -1)   for(j = ii; j < i; j++) sum -= a(i,j) * b(j);
      else if(sum) ii = i;
     b(i) = sum;
   }
   for(i = n - 1; i >= 0; i--) {
     sum = b(i);
     for(j = i+1; j < n; j++) sum -= a(i,j) * b(j);
     b(i) = sum/a(i,i);
   }
} // End: function lubksb()



// Function to invert a matrix
void inverse(Matrix &a, int n)
{        
  int          i,j;
  Vector  indx(n), col(n);
  double       d; 
  Matrix y(n,n);
   
  ludcmp(a, n, indx, &d);   // LU decompose  a[][] 

  // find inverse of a[][] by columns 

  for(j = 0; j < n; j++) {

    // initialize right-side of linear equations 

    for(i = 0; i < n; i++) col(i) = 0.0;
    col(j) = 1.0;

    lubksb(a, n, indx, col);

    // save result in y[][] 

    for(i = 0; i < n; i++) y(i,j) = col(i);
  }   //j-loop over columns 
   
  // return the inverse matrix in a(:,:)

  a = y;
  for(i = 0; i < n; i++) {
    for(j = 0; j < n; j++) a(i,j) = y(i,j);
  } 
}  // End: function inverse()

//Jastrow-faktor
double jastrow(Matrix &R, double beta, double Z, double Z2) {

  double arg = 0;
  for (int i = 1; i < Z; i++)
    for (int j = 0; j < i; j++)
      if ((i < Z2 && j < Z2) || (i >= Z2 && j >= Z2)) {
	double rij = r_ij(R,i,j);
	arg += .25*rij / (1+beta*rij); //samme spinn
      }
      else {
	double rij = r_ij(R,i,j);
	arg += .5*rij / (1+beta*rij); //motsatt spinn
      }
  return exp(arg);
}

//Sjekker om vi f?r singularitetsproblemer med R
bool singularity(Matrix &R, int Z) {

  for (int i = 0; i < Z; i++)
    if (r_i(R,i) < 1e-10)
      return true;

  for (int i = 0; i < Z - 1; i++)
    for (int j = i+1; j < Z; j++)
      if (r_ij(R,i,j) < 1e-10)
	return true;
  return false;
}

//Local energy
double E_local(Matrix &R, double alpha, double beta, int Z, Matrix &F, Matrix & DinvUp,
	       Matrix & DinvDown, int Z2, Matrix &detgrad, Matrix &jastgrad) {

  //Kinetic energy
  double kinetic = 0;
  //Determinantdel
  for (int i = 0; i < Z; i++) {
    for (int j = 0; j < 3; j++) {
      
      if (i < Z2)
	for (int l = 0; l < Z2; l++)
	  kinetic -= phi_deriv2(l,R,i,j,alpha)*DinvUp(l,i);
      else
	for (int l = 0; l < Z2; l++)
	  kinetic -= phi_deriv2(l,R,i,j,alpha)*DinvDown(l,i-Z2);
    }
  }

  //Jastrow part
  double rij,a;
  for (int i = 0; i < Z; i++) {
    for (int j = 0; j < 3; j++) {
      kinetic -= jastgrad(i,j)*jastgrad(i,j);
    }
  }
  for (int i = 0; i < Z-1; i++) {
    for (int j = i+1; j < Z; j++) {
      if ((j < Z2 && i < Z2) || (j >= Z2 && i >= Z2))
	a = .25;
      else
	a = .5;
      rij = r_ij(R,i,j);
      kinetic -= 4*a / (rij*pow(1+beta*rij,3));
    }
  }
  
  //"Interference"-part
  for (int i = 0; i < Z; i++) {
    for (int j = 0; j < 3; j++) {
      kinetic -= 2*detgrad(i,j)*jastgrad(i,j);
    }
  }
  kinetic *= .5;
  //Potential energy, contribution from atomic nucleus
  double potential = 0;
  for (int i = 0; i < Z; i++)
    potential -= Z / r_i(R,i);

  //Electron-electron-potential
  for (int i = 0; i < Z - 1; i++)
    for (int j = i+1; j < Z; j++)
      potential += 1 / r_ij(R,i,j);
  return potential + kinetic;
}

//Function for the quantum force
void calcQF(Matrix &R, Matrix &F, double alpha, double beta, 
	    int Z, Matrix &DinvUp, Matrix &DinvDown, int Z2, Matrix &detgrad, Matrix &jastgrad) {
  
  double sum;
  //Determinant part
  for (int i = 0; i < Z; i++) {
    for (int j = 0; j < 3; j++) {
      sum = 0;
      if (i < Z2)
	for (int l = 0; l < Z2; l++)
	  sum += phi_deriv(l,R,i,j,alpha)*DinvUp(l,i);
      else
	for (int l = 0; l < Z2; l++)
	  sum += phi_deriv(l,R,i,j,alpha)*DinvDown(l,i-Z2);
      detgrad(i,j) = sum;
    }
  }
  //Jastrow part
  double ril,a;
  for (int i = 0; i < Z; i++) {
    for(int j = 0; j < 3; j++) {
      sum = 0;
      for (int l = 0; l < Z; l++) {
	if (l != i) {
	  if ((l < Z2 && i < Z2) || (l >= Z2 && i >= Z2))
	    a = .25;
	  else
	    a = .5;
	  ril = r_ij(R,i,l);
	  sum += (R(i,j)-R(l,j))*a / (ril*pow(1+beta*ril,2));
	}
      }
      jastgrad(i,j) = sum;
    }
  }
  for (int i = 0; i < Z; i++)
    for(int j = 0; j < 3; j++)
      F(i,j) = 2*(detgrad(i,j) + jastgrad(i,j));
}

// Function which implements the Metropolis-Hastings algo and computes the expectation
// values <H> and <H^2>
void MonteCarlo(double &sum, double &sqrsum, int antMC, 
	     double delta_t, long& idum, double alpha, int& N, double beta, int Z,
	     double& sum1alpha, double& sum2alpha, double& sum1beta, double& sum2beta, bool calcderiv) {

  //Initialiser posisjon
  Matrix R(Z,3);
  for (int i = 0; i < Z; i++)
    for (int j = 0; j < 3; j++)
      R(i,j) = gaussian_deviate(&idum);

  int Z2 = Z/2; //dimensjonen til Slater-matrisene

  //Initialize the inverse Slater matrices for spun up and spin down
  Matrix DinvUp(Z2,Z2), DinvDown(Z2,Z2); 
  for (int i = 0; i < Z2; i++) {
    for (int j = 0; j < Z2; j++) {
      DinvUp(i,j) = phi(j,R,i,alpha);
      DinvDown(i,j) = phi(j,R,i+Z2,alpha);
    }
  }
  inverse(DinvUp,Z2);
  inverse(DinvDown,Z2);

  //Inverse Slater matrix for the new position
  Matrix DinvUp_new(Z2,Z2), DinvDown_new(Z2,Z2); 
  for (int i = 0; i < Z2; i++) {
    for (int j = 0; j < Z2; j++) {
      DinvUp_new(i,j) = DinvUp(i,j);
      DinvDown_new(i,j) = DinvDown(i,j);
    }
  }

  //Gradients of the determinants and the Jastrow factor
  Matrix detgrad(Z,3), jastgrad(Z,3), detgrad_new(Z,3), jastgrad_new(Z,3);
  //Initialize the quantum force
  Matrix F(Z,3);
  calcQF(R,F,alpha,beta,Z,DinvUp,DinvDown,Z2,detgrad,jastgrad);
  
  double EL; //Local energy
  double sqrtdt = sqrt(delta_t);
  double D = .5; //The diffusion constant
  double Sj;

  //For Metropolis-Hastings algo:
  Matrix R_new(Z,3), F_new(Z,3);

  double greensratio; // Ratio between Green's functions
  double detratio; //Ratio between Slater determinants in new and old position
  double jastratio; //Ratio between Jastrow factors in new and old position
  double rold,rnew,a;
  double alphaderiv,betaderiv;

  //Loop over Monte Carlo cycles
  for (int k = 0; k < antMC; k++) {

    //We move one electron at the time, this is its loop
    for (int i = 0; i < Z; i++) {

      //move electron i
      for (int j = 0; j < 3; j++)
	R_new(i,j) = R(i,j) + D*delta_t*F(i,j) + gaussian_deviate(&idum)*sqrtdt;
 
      //Leave the other electrons where they are
      for (int ii = 0; ii < Z; ii++)
	for (int j = 0; j < 3; j++)
	  if (ii != i)
	    R_new(ii,j) = R(ii,j);

      //The ratio between Slater determinants
      if (i < Z2) {
	detratio = 0;
	for (int l = 0; l < Z2; l++)
	  detratio += phi(l,R_new,i,alpha) * DinvUp(l,i);
      }
      else {
	detratio = 0;
	for (int l = 0; l < Z2; l++)
	  detratio += phi(l,R_new,i,alpha) * DinvDown(l,i-Z2);
      }

      //Inverse Slater matrix in new position
      if (i < Z2) { //Spin up
	for (int j = 0; j < Z2; j++) {
	  if (j != i) {
	    Sj = 0;
	    for (int l = 0; l < Z2; l++) {
	      Sj += phi(l,R_new,i,alpha) * DinvUp(l,j);
	    }
	    for (int l = 0; l < Z2; l++) 
	      DinvUp_new(l,j) = DinvUp(l,j) - Sj * DinvUp(l,i) / detratio;
	  }
	}
	for (int l = 0; l < Z2; l++)
	  DinvUp_new(l,i) = DinvUp(l,i) / detratio;
      }
      else { //Spin down
	for (int j = 0; j < Z2; j++) {
	  if (j != i-Z2) {
	    Sj = 0;
	    for (int l = 0; l < Z2; l++) {
	      Sj += phi(l,R_new,i,alpha) * DinvDown(l,j);
	    }
	    for (int l = 0; l < Z2; l++) 
	      DinvDown_new(l,j) = DinvDown(l,j) - Sj * DinvDown(l,i-Z2) / detratio;
	  }
	}
	for (int l = 0; l < Z2; l++)
	  DinvDown_new(l,i-Z2) = DinvDown(l,i-Z2) / detratio;
      }

      //Regn ut forhold mellom Jastrow-faktorer
      jastratio = 0;
      for (int l = 0; l < Z; l++) {
	if (l != i) {
	  if ((l < Z2 && i < Z2) || (l >= Z2 && i >= Z2))
	    a = .25;
	  else
	    a = .5;
	  rold = r_ij(R,l,i);
	  rnew = r_ij(R_new,l,i);
	  jastratio += a * (rnew/(1+beta*rnew) - rold/(1+beta*rold));
	}
      }
      jastratio = exp(jastratio);

      //compute the quantum force in the new position
      calcQF(R_new,F_new,alpha,beta,Z,DinvUp_new,DinvDown_new,Z2,detgrad_new,jastgrad_new);

      //compute the ratio between Green's function
      greensratio = 0;
      for (int ii = 0; ii < Z; ii++)
	for (int j = 0; j < 3; j++)
	  greensratio += .5*(F_new(ii,j)+F(ii,j)) * (.5*D*delta_t*(F(ii,j)-F_new(ii,j)) + R(ii,j) - R_new(ii,j));
      greensratio = exp(greensratio);
	
      
      //Metropolis-Hastings test
      if (ran1(&idum) < greensratio*detratio*detratio*jastratio*jastratio) {
	//Accept the new position:
	//Update the Slater matrix
	if (i < Z2)
	  for (int l = 0; l < Z2; l++)
	    for (int m = 0; m < Z2; m++)
	      DinvUp(l,m) = DinvUp_new(l,m);
	else
	  for (int l = 0; l < Z2; l++)
	    for (int m = 0; m < Z2; m++)
	      DinvDown(l,m) = DinvDown_new(l,m);

	//Update positions, quantum force and gradients
	for (int ii = 0; ii < Z; ii++) {
	  for (int j = 0; j < 3; j++) {
	    R(ii,j) = R_new(ii,j);
	    F(ii,j) = F_new(ii,j);
	    detgrad(ii,j) = detgrad_new(ii,j);
	    jastgrad(ii,j) = jastgrad_new(ii,j);
	  }
	}

      }

    }
    
    if (!singularity(R,Z) && k > antMC/20) { 

      EL = E_local(R,alpha,beta,Z,F,DinvUp,DinvDown,Z2,detgrad,jastgrad);
      
      sum += EL;
      sqrsum += EL * EL;
      if (calcderiv) {
	alphaderiv = (slater(R,alpha+h,Z,Z2) - slater(R,alpha-h,Z,Z2))/(2*h*slater(R,alpha,Z,Z2));
	sum1alpha += alphaderiv*EL;
	sum2alpha += alphaderiv;
	betaderiv = (jastrow(R,beta+h,Z,Z2) - jastrow(R,beta-h,Z,Z2))/(2*h*jastrow(R,beta,Z,Z2));
	sum1beta += betaderiv*EL;
	sum2beta += betaderiv;
      }
      N++;
    }
  } //End loop over Monte Carlo cycles

}

//Energifunksjon for dfpmin
double energy(Vector &param) {

  double sum = 0;
  double sqrsum = 0;
  int N = 0;
  double sum1alpha = 0;
  double sum2alpha = 0;
  double sum1beta = 0;
  double sum2beta = 0;
  if (param(1) < 0 && Z_global == 4)
    return -12;
  if (param(1) < 0 && Z_global == 10)
    return -100;
  MonteCarlo(sum,sqrsum,antMC_global,delta_t_global,idum_global,param(0),N,param(1),
	  Z_global,sum1alpha,sum2alpha,sum1beta,sum2beta,false);
  cout << "alpha = " << param(0) << "   beta =  " << param(1) << "    energy = " << sum/N << endl;

  return sum / N;
}

//Derivertfunksjon for dfpmin
void energy_deriv(Vector &param, Vector &values) {
  
  double sum = 0;
  double sqrsum = 0;
  int N = 0;
  double sum1alpha = 0;
  double sum2alpha = 0;
  double sum1beta = 0;
  double sum2beta = 0;

  if (param(0) < 0 && Z_global == 4) {
    values(0) = .1;
    values(1) = -.5;
    return;
  }
  if (param(1) < 0 && Z_global == 10) {
    values(0) = .1;
    values(1) = -10;
    return;
  }

  MonteCarlo(sum,sqrsum,antMC_global,delta_t_global,idum_global,param(0),N,param(1),
	  Z_global,sum1alpha,sum2alpha,sum1beta,sum2beta,true);
  
  double E = sum / N;
  
  values(0) = 2 * ( (sum1alpha/N) - (sum2alpha/N)*E );
  values(1) = 2 * ( (sum1beta/N) - (sum2beta/N)*E );

  cout << "alpha = " << param(0) << "   beta =  " << param(1) << 
    "    alphaderiv = " << values(0) << "    betaderiv = " << values(1) << endl;
}

int main() {

  //MPI-initialisering
  int numprocs;
  int minrank;
  //  MPI_Init(NULL,NULL);
  //MPI_Comm_size(MPI_COMM_WORLD, &numprocs);
  //MPI_Comm_rank(MPI_COMM_WORLD, &minrank);


  int antMC = 5000000; //number of cycles
  antMC_global = antMC;
  int Z = 4; // Charge of atomic nucleus
  Z_global = Z;
  double alpha = 3.6 + minrank*0.2; 
  double beta = .1 + minrank*0.1; 
  long idum = - time(NULL) +  minrank*1000;
  idum_global = idum;
  double delta_t = .05; //Time step
  delta_t_global = delta_t;

  Vector param(2), values(2);
  param(0) = alpha;
  param(1) = beta;

  double gtol = .001;
  int iter;
  double fret;

  dfpmin(param, 2, gtol, &iter, &fret, &energy, &energy_deriv);

  cout << "alpha = " << param(0) << endl;
  cout << "beta = " << param(1) << endl;
  cout << "energy = " << fret << endl;
  cout << "iter = " << iter << endl << endl;
    
  //Avslutt MPI
  //  MPI_Finalize();

  return 0;
}