C++ Minimizer类代码示例

double run_soft_sphere(double reduced_density, double temp) {
  Functional f = SoftFluid(sigma, 1, 0);
  const double mu = find_chemical_potential(OfEffectivePotential(f), temp, reduced_density*pow(2,-5.0/2.0));
  printf("mu is %g for reduced_density = %g at temperature %g\n", mu, reduced_density, temp);

  //printf("Filling fraction is %g with functional %s at temperature %g\n", reduced_density, teff);
  //fflush(stdout);
  temperature = temp;
  //if (kT == 0) kT = ;1

  Lattice lat(Cartesian(xmax,0,0), Cartesian(0,ymax,0), Cartesian(0,0,zmax));
  GridDescription gd(lat, dx);

  Grid softspherepotential(gd);
  softspherepotential.Set(soft_sphere_potential);

  f = SoftFluid(sigma, 1, mu); // compute approximate energy with chemical potential mu
  const double approx_energy = f(temperature, reduced_density*pow(2,-5.0/2.0))*xmax*ymax*zmax;
  const double precision = fabs(approx_energy*1e-9);

  f = OfEffectivePotential(SoftFluid(sigma, 1, mu) + ExternalPotential(softspherepotential));

  static Grid *potential = 0;
  potential = new Grid(gd);
  *potential = softspherepotential - temperature*log(reduced_density*pow(2,-5.0/2.0)/(1.0*radius*radius*radius))*VectorXd::Ones(gd.NxNyNz); // Bad starting guess
  printf("\tMinimizing to %g absolute precision from %g from %g...\n", precision, approx_energy, temperature);
  fflush(stdout);

  Minimizer min = Precision(precision,
                            PreconditionedConjugateGradient(f, gd, temperature,
                                potential,
                                QuadraticLineMinimizer));
  took("Setting up the variables");
  for (int i=0; min.improve_energy(true) && i<100; i++) {
  }

  took("Doing the minimization");
  min.print_info();

  Grid density(gd, EffectivePotentialToDensity()(temperature, gd, *potential));
  //printf("# per area is %g at filling fraction %g\n", density.sum()*gd.dvolume/dw/dw, reduced_density);

  char *plotname = (char *)malloc(1024);

  sprintf(plotname, "papers/fuzzy-fmt/figs/radial-wca-%06.4f-%04.2f.dat", temp, reduced_density);
  z_plot(plotname, Grid(gd, pow(2,5.0/2.0)*density));
  free(plotname);

  {
    //double peak = peak_memory()/1024.0/1024;
    //double current = current_memory()/1024.0/1024;
    //printf("Peak memory use is %g M (current is %g M)\n", peak, current);

  }

  took("Plotting stuff");
  printf("density %g gives ff %g for reduced_density = %g and T = %g\n", density(0,0,gd.Nz/2),
         density(0,0,gd.Nz/2)*4*M_PI/3, reduced_density, temp);
  return density(0, 0, gd.Nz/2)*4*M_PI/3; // return bulk filling fraction
}

int main(int, char **) {
  FILE *o = fopen("paper/figs/constrained-water.dat", "w");

  Functional f = OfEffectivePotential(SaftFluidSlow(water_prop.lengthscale,
                                                    water_prop.epsilonAB, water_prop.kappaAB,
                                                    water_prop.epsilon_dispersion,
                                                    water_prop.lambda_dispersion, water_prop.length_scaling, 0));
  double mu_satp = find_chemical_potential(f, water_prop.kT,
                                           water_prop.liquid_density);
  Lattice lat(Cartesian(width,0,0), Cartesian(0,width,0), Cartesian(0,0,zmax));
  GridDescription gd(lat, 0.1);

  Grid potential(gd);
  Grid constraint(gd);
  constraint.Set(notinwall);

  f = constrain(constraint,
                OfEffectivePotential(SaftFluidSlow(water_prop.lengthscale,
                                                   water_prop.epsilonAB, water_prop.kappaAB,
                                                   water_prop.epsilon_dispersion,
                                                   water_prop.lambda_dispersion, water_prop.length_scaling, mu_satp)));


  Minimizer min = Precision(0, PreconditionedConjugateGradient(f, gd, water_prop.kT, &potential,
                                                               QuadraticLineMinimizer));

  potential = water_prop.liquid_density*constraint
    + water_prop.vapor_density*VectorXd::Ones(gd.NxNyNz);
  //potential = water_prop.liquid_density*VectorXd::Ones(gd.NxNyNz);
  potential = -water_prop.kT*potential.cwise().log();

  const int numiters = 50;
  for (int i=0;i<numiters && min.improve_energy(true);i++) {
    fflush(stdout);
    Grid density(gd, EffectivePotentialToDensity()(water_prop.kT, gd, potential));
    density.epsNative1d("paper/figs/1d-constrained-plot.eps",
			Cartesian(0,0,0), Cartesian(0,0,zmax),
			water_prop.liquid_density, 1, "Y axis: , x axis: ");
  }
  min.print_info();

  double energy = min.energy()/width/width;
  printf("Energy is %.15g\n", energy);

  double N = 0;
  {
    Grid density(gd, EffectivePotentialToDensity()(water_prop.kT, gd, potential));
    for (int i=0;i<gd.NxNyNz;i++) N += density[i]*gd.dvolume;
  }
  N = N/width/width;
  printf("N is %.15g\n", N);

    Grid density(gd, EffectivePotentialToDensity()(water_prop.kT, gd, potential));
    density.epsNative1d("paper/figs/1d-constrained-plot.eps", Cartesian(0,0,0), Cartesian(0,0,zmax), water_prop.liquid_density, 1, "Y axis: , x axis: ");
    //potential.epsNative1d("hard-wall-potential.eps", Cartesian(0,0,0), Cartesian(0,0,zmax), 1, 1);

    fclose(o);
}

double run_walls(double reduced_density, const char *name, Functional fhs, double teff) {
  double kT = teff;
  if (kT == 0) kT = 1;

  Functional f = OfEffectivePotential(fhs);

  const double zmax = width + 2*spacing;
  Lattice lat(Cartesian(dw,0,0), Cartesian(0,dw,0), Cartesian(0,0,zmax));
  GridDescription gd(lat, dx);

  Grid constraint(gd);
  constraint.Set(notinwall);
  f = constrain(constraint, f);

  Grid potential(gd);
  potential = pow(2,-5.0/2.0)*(reduced_density*constraint + 1e-4*reduced_density*VectorXd::Ones(gd.NxNyNz));
  potential = -kT*potential.cwise().log();

  const double approx_energy = fhs(kT, reduced_density*pow(2,-5.0/2.0))*dw*dw*width;
  const double precision = fabs(approx_energy*1e-11);
  printf("\tMinimizing to %g absolute precision from %g from %g...\n", precision, approx_energy, kT);
  fflush(stdout);

  Minimizer min = Precision(precision,
                            PreconditionedConjugateGradient(f, gd, kT,
                                                            &potential,
                                                            QuadraticLineMinimizer));
  took("Setting up the variables");
  if (strcmp(name, "hard") != 0 && false) {
    printf("For now, SoftFluid doesn't work properly, so we're skipping the\n");
    printf("minimization at temperature %g.\n", teff);
  } else {
    for (int i=0;min.improve_energy(false) && i<100;i++) {
    }
  }
  took("Doing the minimization");
  min.print_info();

  Grid density(gd, EffectivePotentialToDensity()(kT, gd, potential));
  //printf("# per area is %g at filling fraction %g\n", density.sum()*gd.dvolume/dw/dw, eta);

  char *plotname = (char *)malloc(1024);

  sprintf(plotname, "papers/fuzzy-fmt/figs/walls%s-%06.4f-%04.2f.dat", name, teff, reduced_density);
  z_plot(plotname, Grid(gd, density*pow(2,5.0/2.0)));
  free(plotname);

  took("Plotting stuff");
  printf("density %g gives ff %g for reduced density = %g and T = %g\n", density(0,0,gd.Nz/2),
         density(0,0,gd.Nz/2)*4*M_PI/3, reduced_density, teff);
  return density(0, 0, gd.Nz/2)*4*M_PI/3; // return bulk filling fraction
}

void TForm2::testMinimize() {
     		double *params = new1D(2, 0.0);
		MinimizerFunction *mf = new MinimizerFunction(params, 2);
		Minimizer *m = new Minimizer(mf);
		m->DoMinimize();
		double* pOut = mf->getParameters();

		p1Edit->Text = UnicodeString(pOut[0]);
		p2Edit->Text = UnicodeString(pOut[1]);

		delete m;
		delete mf;
		delete1D(params);
}

bool ConjugateGradientType::improve_energy(bool verbose) {
  iter++;
  //printf("I am running ConjugateGradient::improve_energy\n");
  const double E0 = energy();
  if (E0 != E0) {
    // There is no point continuing, since we're starting with a NaN!
    // So we may as well quit here.
    if (verbose) {
      printf("The initial energy is a NaN, so I'm quitting early from ConjugateGradientType::improve_energy.\n");
      f.print_summary("has nan:", E0);
      fflush(stdout);
    }
    return false;
  }
  double gdotd;
  {
    const VectorXd g = -grad();
    // Let's immediately free the cached gradient stored internally!
    invalidate_cache();

    // Note: my notation vaguely follows that of
    // [wikipedia](http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method).
    // I use the Polak-Ribiere method, with automatic direction reset.
    // Note that we could save some memory by using Fletcher-Reeves, and
    // it seems worth implementing that as an option for
    // memory-constrained problems (then we wouldn't need to store oldgrad).
    double beta = g.dot(g - oldgrad)/oldgradsqr;
    oldgrad = g;
    if (beta < 0 || beta != beta || oldgradsqr == 0) beta = 0;
    oldgradsqr = oldgrad.dot(oldgrad);
    direction = g + beta*direction;
    gdotd = oldgrad.dot(direction);
    if (gdotd < 0) {
      direction = oldgrad; // If our direction is uphill, reset to gradient.
      if (verbose) printf("reset to gradient...\n");
      gdotd = oldgrad.dot(direction);
    }
  }

  Minimizer lm = linmin(f, gd, kT, x, direction, -gdotd, &step);
  for (int i=0; i<100 && lm.improve_energy(verbose); i++) {
    if (verbose) lm.print_info("\t");
  }
  if (verbose) {
    //lm->print_info();
    print_info();
    printf("grad*dir/oldgrad*dir = %g\n", grad().dot(direction)/gdotd);
  }
  return (energy() < E0);
}

bool PreconditionedConjugateGradientType::improve_energy(bool verbose) {
  iter++;
  //printf("I am running ConjugateGradient::improve_energy\n");
  const double E0 = energy();
  if (isnan(E0)) {
    // There is no point continuing, since we're starting with a NaN!
    // So we may as well quit here.
    if (verbose) {
      printf("The initial energy is a NaN, so I'm quitting early.\n");
      fflush(stdout);
    }
    return false;
  }
  double beta;
  {
    // Note: my notation vaguely follows that of
    // [wikipedia](http://en.wikipedia.org/wiki/Nonlinear_conjugate_gradient_method).
    // I use the Polak-Ribiere method, with automatic direction reset.
    // Note that we could save some memory by using Fletcher-Reeves, and
    // it seems worth implementing that as an option for
    // memory-constrained problems (then we wouldn't need to store oldgrad).
    pgrad(); // compute pgrad first, since that computes both.
    beta = -pgrad().dot(-grad() - oldgrad)/oldgradsqr;
    oldgrad = -grad();
    if (beta < 0 || beta != beta || oldgradsqr == 0) beta = 0;
    if (verbose) printf("beta = %g\n", beta);
    oldgradsqr = -pgrad().dot(oldgrad);
    direction = -pgrad() + beta*direction;
    // Let's immediately free the cached gradient stored internally!
    invalidate_cache();
  } // free g and pg!

  const double gdotd = oldgrad.dot(direction);

  Minimizer lm = linmin(f, gd, kT, x, direction, -gdotd, &step);
  for (int i=0; i<100 && lm.improve_energy(verbose); i++) {
    if (verbose) lm.print_info("\t");
  }
  if (verbose) {
    //lm->print_info();
    print_info();
    printf("grad*oldgrad = %g\n", grad().dot(direction)/gdotd);
  }
  return (energy() < E0 || beta != 0);
}

int main (int argc, char ** argv) {
	H.init ("Test: circle packing", argc, argv, "s=4");

	Map m (13);
	m	<< Edge(0,1) << Edge(0,3) << Edge(0,5)
	 	<< Edge(1,2) << Edge(1,4) << Edge(1,6) << Edge(1,3) << Edge(1,0)
	 	<< Edge(2,7) << Edge(2,4) << Edge(2,1)
	 	<< Edge(3,0) << Edge(3,1) << Edge(3,6) << Edge(3,5)
	 	<< Edge(4,1) << Edge(4,2) << Edge(4,7) << Edge(4,6)
	 	<< Edge(5,0) << Edge(5,3) << Edge(5,6) << Edge(5,8) << Edge(5,10)
	 	<< Edge(6,1) << Edge(6,4) << Edge(6,7) << Edge(6,9) << Edge(6,11) << Edge(6,8) << Edge(6,5) << Edge(6,3)
	 	<< Edge(7,12) << Edge(7,9) << Edge(7,6) << Edge(7,4) << Edge(7,2)
	 	<< Edge(8,5) << Edge(8,6) << Edge(8,11) << Edge(8,10)
	 	<< Edge(9,6) << Edge(9,7) << Edge(9,12) << Edge(9,11)
	 	<< Edge(10,5) << Edge(10,8) << Edge(10,11)
	 	<< Edge(11,10) << Edge(11,8) << Edge(11,6) << Edge(11,9) << Edge(11,12)
	 	<< Edge(12,11) << Edge(12,9) << Edge(12,7);

	for (int i=0; i<int(H['s']); ++i) m.barycentric();
	m.inscribe(m.face(Edge(1,m.v[1]->adj.back()))); m.balance(); m.show();

	Vector<double> x(3*m.n); double r = 1.0/sqrt(m.n);

	for (int i=0; i<m.n; ++i) {
		x[3*i]        = m.v[i]->z.real() / (1-r);
		x[3*i+1]      = m.v[i]->z.imag() / (1-r);
		x[3*i+2]      = .8*r;
	}

	Minimizer<double> MM (3*m.n, [&m](const Vector<double> &x, Vector<double> &g) { return Map_fg_circle_disk (x,g,&m); }); MM.cb = cb;
	MM.minimize_qn (x); x = MM.x;

	std::cerr << "Number of vertices:    " << m.n << std::endl;
	std::cerr << "Final value of f:      " << MM.fx << std::endl;
	std::cerr << "Final square gradient: " << inner_prod(MM.gx,MM.gx) << std::endl;

	for (int i=0; i<m.n; ++i) {
		m.v[i]->z = cpx (x[3*i], x[3*i+1]);
		m.v[i]->r = x[3*i+2];
	}

	Figure f; m.plot_circles(f); f.add (new Circle(cpx(0.0,0.0),1.0)); f.show(); f.pause(); f.output ();
}

double surface_tension(Minimizer min, Functional f0, LiquidProperties prop,
                       bool verbose, const char *plotname) {
  int numptspersize = 100;
  int size = 64;
  const int gas_size = 10;
  Lattice lat(Cartesian(1,0,0), Cartesian(0,1,0), Cartesian(0,0,size*prop.lengthscale));
  GridDescription gd(lat, 1, 1, numptspersize*size);
  Grid potential(gd);

  // Set the density to range from vapor to liquid
  const double Veff_liquid = -prop.kT*log(prop.liquid_density);
  const double Veff_gas = -prop.kT*log(prop.vapor_density);
  for (int i=0; i<gd.NxNyNz*gas_size/size; i++) potential[i] = Veff_gas;
  for (int i=gd.NxNyNz*gas_size/size; i<gd.NxNyNz; i++) potential[i] = Veff_liquid;

  // Enable the following line for debugging...
  //f0.run_finite_difference_test("f0", prop.kT, potential);
  min.minimize(f0, gd, &potential);
  while (min.improve_energy(verbose))
    if (verbose) {
      printf("Working on liberated interface...\n");
      fflush(stdout);
    }
  const double Einterface = f0.integral(prop.kT, potential);
  double Ninterface = 0;
  Grid interface_density(gd, EffectivePotentialToDensity()(prop.kT, gd, potential));
  for (int i=0;i<gd.NxNyNz;i++) Ninterface += interface_density[i]*gd.dvolume;
  if (verbose) printf("Got interface energy of %g.\n", Einterface);
  
  for (int i=0; i<gd.NxNyNz; i++) potential[i] = Veff_gas;
  min.minimize(f0, gd, &potential);
  while (min.improve_energy(verbose))
    if (verbose) {
      printf("Working on gas...\n");
      fflush(stdout);
    }
  const double Egas = f0.integral(prop.kT, potential);
  double Ngas = 0;
  {
    Grid density(gd, EffectivePotentialToDensity()(prop.kT, gd, potential));
    for (int i=0;i<gd.NxNyNz;i++) Ngas += density[i]*gd.dvolume;
  }
  
  for (int i=0; i<gd.NxNyNz; i++) potential[i] = Veff_liquid;
  if (verbose) {
    printf("\n\n\nWorking on liquid...\n");
    fflush(stdout);
  }
  min.minimize(f0, gd, &potential);
  while (min.improve_energy(verbose))
    if (verbose) {
      printf("Working on liquid...\n");
      fflush(stdout);
    }
  const double Eliquid = f0.integral(prop.kT, potential);
  double Nliquid = 0;
  {
    Grid density(gd, EffectivePotentialToDensity()(prop.kT, gd, potential));
    for (int i=0;i<gd.NxNyNz;i++) Nliquid += density[i]*gd.dvolume;
  }
  
  const double X = Ninterface/Nliquid; // Fraction of volume effectively filled with liquid.
  const double surface_tension = (Einterface - Eliquid*X - Egas*(1-X))/2;
  if (verbose) {
    printf("\n\n");
    printf("interface energy is %.15g\n", Einterface);
    printf("gas energy is %.15g\n", Egas);
    printf("liquid energy is %.15g\n", Eliquid);
    printf("Ninterface/liquid/gas = %g/%g/%g\n", Ninterface, Nliquid, Ngas);
    printf("X is %g\n", X);
    printf("surface tension is %.10g\n", surface_tension);
  }
  if (plotname)
    interface_density.Dump1D(plotname, Cartesian(0,0,0),
                             Cartesian(0,0,size*prop.lengthscale));
  return surface_tension;
}

//.........这里部分代码省略.........
    } else if (run == "N2"){

        double d = 2.116115162;
        rowvec posN1 = {-0.5*d, 0.0, 0.0};
        rowvec posN2= {0.5*d, 0.0, 0.0};
        rowvec charges = {7.0, 7.0};
        int nElectrons = 14;

        mat nucleiPositions = zeros<mat>(2,3);
        nucleiPositions.row(0) = posN1;
        nucleiPositions.row(1) = posN2;

        BasisFunctions* basisFunctions = new BasisFunctions;
        basisFunctions->addContracteds("../../HartreeFock/inFiles/basisSets/N_6311Gs.dat", 0);
        basisFunctions->addContracteds("../../HartreeFock/inFiles/basisSets/N_6311Gs.dat", 1);

        System *system;
        system = new System(basisFunctions, nucleiPositions, charges, nElectrons);

        RMP solver(system,2);
        solver.solve();

        if (my_rank == 0){
            cout << "Energy: " << setprecision(9) << solver.getEnergyHF() + solver.getEnergy2order() + solver.getEnergy3order() << endl;
        }

        delete system;
        delete basisFunctions;

    } else if (run == "FCl") {

        double d = 3.154519593;
        rowvec posF = {-0.5*d, 0.0, 0.0};
        rowvec posCl= {0.5*d, 0.0, 0.0};
        rowvec charges = {9.0, 17.0};
        int nElectrons = 26;

        mat nucleiPositions = zeros<mat>(2,3);
        nucleiPositions.row(0) = posF;
        nucleiPositions.row(1) = posCl;

        BasisFunctions* basisFunctions = new BasisFunctions;
        basisFunctions->addContracteds("../../HartreeFock/inFiles/basisSets/F_6311Gs.dat", 0);
        basisFunctions->addContracteds("../../HartreeFock/inFiles/basisSets/Cl_6311Gs.dat", 1);

        System *system;
        system = new System(basisFunctions, nucleiPositions, charges, nElectrons);

        RMP solver(system,2);
        solver.solve();

        if (my_rank == 0){
            cout << "Energy: " << setprecision(9) << solver.getEnergyHF() + solver.getEnergy2order() << endl;
        }

    } else if (run == "H2O_Minimize"){

        rowvec O = {0.0,0.0,0.0};
        rowvec H1 = {1.0,0.0,0.0};
        rowvec H2 = {0.0,1.0,0.0};
        rowvec charges = {8.0,1.0,1.0};
        int nElectrons = 10;
        mat nucleiPositions = zeros<mat>(3,3);
        nucleiPositions.row(0) = O;
        nucleiPositions.row(1) = H1;
        nucleiPositions.row(2) = H2;
        BasisFunctions* basisFunctions = new BasisFunctions;
        basisFunctions->addContracteds("../../HartreeFock/inFiles/basisSets/O_431G.dat", 0);
        basisFunctions->addContracteds("../../HartreeFock/inFiles/basisSets/H_431G.dat", 1);
        basisFunctions->addContracteds("../../HartreeFock/inFiles/basisSets/H_431G.dat", 2);
        System *system = new System(basisFunctions, nucleiPositions, charges, nElectrons);
        RMP *solver = new RMP(system,1);
        HartreeFockFunc *func = new HartreeFockFunc(solver, system);
        Minimizer *minimizer = new Minimizer(func);
        minimizer->solve();

        if (my_rank == 0){
            cout << system->getNucleiPositions() << endl;
            cout << setprecision(7) << system->getNucleiPositions()(1,0) << endl;
            cout << "Energy min.: " << setprecision(14) << minimizer->getMinValue() << endl;
            cout << "Energy max.: " << setprecision(14) << minimizer->getMaxValue() << endl;
        }

    } else {
        if (my_rank == 0){
            cout << "No valid run selected." << endl;
        }
    }

    clock_t end = clock();
    if (my_rank == 0){
        cout << "Elapsed time: "<< (double(end - begin))/CLOCKS_PER_SEC << endl;
    }

#ifdef RUN_MPI
    MPI_Finalize();
#endif

    return 0;
}

int main(int argc, char *argv[]) {
  clock_t start_time = clock();
  if (argc > 1) {
    if (sscanf(argv[1], "%lg", &diameter) != 1) {
      printf("Got bad argument: %s\n", argv[1]);
      return 1;
    }
    diameter *= nm;
    using_default_diameter = false;
  }
  printf("Diameter is %g bohr = %g nm\n", diameter, diameter/nm);
  const double padding = 1*nm;
  xmax = ymax = zmax = diameter + 2*padding;

  char *datname = (char *)malloc(1024);
  sprintf(datname, "papers/hughes-saft/figs/sphere-%04.2fnm-energy.dat", diameter/nm);
  
  Functional f = OfEffectivePotential(SaftFluid2(hughes_water_prop.lengthscale,
						hughes_water_prop.epsilonAB, hughes_water_prop.kappaAB,
						hughes_water_prop.epsilon_dispersion,
						hughes_water_prop.lambda_dispersion,
						hughes_water_prop.length_scaling, 0));
  double n_1atm = pressure_to_density(f, hughes_water_prop.kT, atmospheric_pressure,
				      0.001, 0.01);

  double mu_satp = find_chemical_potential(f, hughes_water_prop.kT, n_1atm);

  f = OfEffectivePotential(SaftFluid2(hughes_water_prop.lengthscale,
				     hughes_water_prop.epsilonAB, hughes_water_prop.kappaAB,
				     hughes_water_prop.epsilon_dispersion,
				     hughes_water_prop.lambda_dispersion,
				     hughes_water_prop.length_scaling, mu_satp));
  Functional S = OfEffectivePotential(EntropySaftFluid2(new_water_prop.lengthscale,
                                                        new_water_prop.epsilonAB,
                                                        new_water_prop.kappaAB,
                                                        new_water_prop.epsilon_dispersion,
                                                        new_water_prop.lambda_dispersion,
                                                        new_water_prop.length_scaling));
  
  const double EperVolume = f(hughes_water_prop.kT, -hughes_water_prop.kT*log(n_1atm));
  const double EperNumber = EperVolume/n_1atm;
  const double SperNumber = S(hughes_water_prop.kT, -hughes_water_prop.kT*log(n_1atm))/n_1atm;
  const double EperCell = EperVolume*(zmax*ymax*xmax - (M_PI/6)*diameter*diameter*diameter);

  //for (diameter=0*nm; diameter<3.0*nm; diameter+= .1*nm) {
    Lattice lat(Cartesian(xmax,0,0), Cartesian(0,ymax,0), Cartesian(0,0,zmax));
    GridDescription gd(lat, 0.2);
    
    Grid potential(gd);
    Grid constraint(gd);
    constraint.Set(notinwall);
    
    f = OfEffectivePotential(SaftFluid2(hughes_water_prop.lengthscale,
				       hughes_water_prop.epsilonAB, hughes_water_prop.kappaAB,
				       hughes_water_prop.epsilon_dispersion,
				       hughes_water_prop.lambda_dispersion,
				       hughes_water_prop.length_scaling, mu_satp));
    f = constrain(constraint, f);
    //constraint.epsNativeSlice("papers/hughes-saft/figs/sphere-constraint.eps",
    // 			      Cartesian(0,ymax,0), Cartesian(0,0,zmax), 
    // 			      Cartesian(0,ymax/2,zmax/2));
    //printf("Constraint has become a graph!\n");
   
    potential = hughes_water_prop.liquid_density*constraint
      + 100*hughes_water_prop.vapor_density*VectorXd::Ones(gd.NxNyNz);
    //potential = hughes_water_prop.liquid_density*VectorXd::Ones(gd.NxNyNz);
    potential = -hughes_water_prop.kT*potential.cwise().log();
    
    double energy;
    {
      const double surface_tension = 5e-5; // crude guess from memory...
      const double surfprecision = 1e-4*M_PI*diameter*diameter*surface_tension; // four digits precision
      const double bulkprecision = 1e-12*fabs(EperCell); // but there's a limit on our precision for small spheres
      const double precision = bulkprecision + surfprecision;
      Minimizer min = Precision(precision,
                                PreconditionedConjugateGradient(f, gd, hughes_water_prop.kT, 
                                                                &potential,
                                                                QuadraticLineMinimizer));
      
      printf("\nDiameter of sphere = %g bohr (%g nm)\n", diameter, diameter/nm);
      
      const int numiters = 200;
      for (int i=0;i<numiters && min.improve_energy(true);i++) {
        //fflush(stdout);
        //Grid density(gd, EffectivePotentialToDensity()(hughes_water_prop.kT, gd, potential));
        
        //density.epsNativeSlice("papers/hughes-saft/figs/sphere.eps", 
        //			     Cartesian(0,ymax,0), Cartesian(0,0,zmax), 
        //			     Cartesian(0,ymax/2,zmax/2));
        
        //sleep(3);

        double peak = peak_memory()/1024.0/1024;
        double current = current_memory()/1024.0/1024;
        printf("Peak memory use is %g M (current is %g M)\n", peak, current);
      }
      
      double peak = peak_memory()/1024.0/1024;
      double current = current_memory()/1024.0/1024;
      printf("Peak memory use is %g M (current is %g M)\n", peak, current);
//.........这里部分代码省略.........

int main(int, char **) {
    for (double eta = 0.3; eta < 0.35; eta += 0.1) {
        // Generates a data file for the pair distribution function, for filling fraction eta
        // and distance of first sphere from wall of z0. Data saved in a table such that the
        // columns are x values and rows are z1 values.
        printf("Now starting sphere_with_wall with eta = %g\n",eta);
        Lattice lat(Cartesian(width,0,0), Cartesian(0,width,0), Cartesian(0,0,width+2*spacing));
        GridDescription gd(lat, dx); //the resolution here dramatically affects our memory use

        Functional f = OfEffectivePotential(WB + IdealGas());
        double mu = find_chemical_potential(f, 1, eta/(4*M_PI/3));
        f = OfEffectivePotential(WB + IdealGas()
                                 + ChemicalPotential(mu));

        Grid potential(gd);
        Grid constraint(gd);

        constraint.Set(*notinwall_or_sphere);
        constraint.epsNativeSlice("myconstraint.eps",
                                  Cartesian(0, 0, 2*(width+2*spacing)),
                                  Cartesian(2*width, 0, 0),
                                  Cartesian(0, 0, 0));

        f = constrain(constraint, f);

        potential = (eta*constraint + 1e-4*eta*VectorXd::Ones(gd.NxNyNz))/(4*M_PI/3);
        potential = -potential.cwise().log();

        const double approx_energy = (WB + IdealGas() + ChemicalPotential(mu))(1, eta/(4*M_PI/3))*dw*dw*width;
        const double precision = fabs(approx_energy*1e-4);
        //printf("Minimizing to %g absolute precision...\n", precision);
        Minimizer min = Precision(precision,
                                  PreconditionedConjugateGradient(f, gd, 1,
                                          &potential,
                                          QuadraticLineMinimizer));
        {
            double peak = peak_memory()/1024.0/1024;
            double current = current_memory()/1024.0/1024;
            printf("Peak memory use is %g M (current is %g M)\n", peak, current);
            fflush(stdout);
        }

        for (int i=0; min.improve_energy(true) && i<100; i++) {
            double peak = peak_memory()/1024.0/1024;
            double current = current_memory()/1024.0/1024;
            printf("Peak memory use is %g M (current is %g M)\n", peak, current);
            fflush(stdout);
        }
        Grid density(gd, EffectivePotentialToDensity()(1, gd, potential));

        char *plotname = new char[1024];
        sprintf(plotname, "papers/pair-correlation/figs/walls/wallsWB-sphere-dft-%04.2f.dat", eta);
        pair_plot(plotname, density);
        delete[] plotname;

        char *plotname_path = new char[1024];
        sprintf(plotname_path, "papers/pair-correlation/figs/walls/wallsWB-sphere-dft-path-%04.2f.dat", eta);
        path_plot(plotname_path, density, constraint);
        delete[] plotname_path;
        fflush(stdout);
        {
            double peak = peak_memory()/1024.0/1024;
            double current = current_memory()/1024.0/1024;
            printf("Peak memory use is %g M (current is %g M)\n", peak, current);
            fflush(stdout);
        }
        fflush(stdout);
    }
    fflush(stdout);
    // Just create this file so make knows we have run.
    if (!fopen("papers/pair-correlation/figs/walls_sphere.dat", "w")) {
        printf("Error creating walls.dat!\n");
        return 1;
    }
    fflush(stdout);
    return 1;
}

int main(int argc, char *argv[]) {
  if (argc == 5) {
    if (sscanf(argv[1], "%lg", &xmax) != 1) {
      printf("Got bad x argument: %s\n", argv[1]);
      return 1;
    }
    if (sscanf(argv[2], "%lg", &ymax) != 1) {
      printf("Got bad y argument: %s\n", argv[2]);
      return 1;
    }
    if (sscanf(argv[3], "%lg", &zmax) != 1) {
      printf("Got bad z argument: %s\n", argv[3]);
      return 1;
    }
    if (sscanf(argv[4], "%lg", &N) != 1) {
      printf("Got bad N argument: %s\n", argv[4]);
      return 1;
    }
    using_default_box = false;
    printf("Box is %g x %g x %g hard sphere diameters, and it holds %g of them\n", xmax, ymax, zmax, N);
  }

  char *datname = (char *)malloc(1024);
  sprintf(datname, "papers/contact/figs/box-%02.0f,%02.0f,%02.0f-%02.0f-energy.dat", xmax, ymax, zmax, N);
  
  FILE *o = fopen(datname, "w");

  const double myvolume = (xmax+2)*(ymax+2)*(zmax+2);
  const double meandensity = N/myvolume;

  Functional f = OfEffectivePotential(HS + IdealGas());
  double mu = find_chemical_potential(f, 1, meandensity);
  f = OfEffectivePotential(HS + IdealGas()
                           + ChemicalPotential(mu));

  Lattice lat(Cartesian(xmax+3,0,0), Cartesian(0,ymax+3,0), Cartesian(0,0,zmax+3));
  GridDescription gd(lat, 0.05);
    
  Grid potential(gd);
  Grid constraint(gd);
  constraint.Set(notinwall);
  took("Setting the constraint");

  printf("xmax = %g\nymax = %g\nzmax = %g\nmeandensity=%g\n", xmax, ymax, zmax, meandensity);
  f = constrain(constraint, f);
  constraint.epsNativeSlice("papers/contact/figs/box-constraint.eps",
   			      Cartesian(0,ymax+4,0), Cartesian(0,0,zmax+4), 
   			      Cartesian(0,-ymax/2-2,-zmax/2-2));
  printf("Constraint has become a graph!\n");
  
  potential = meandensity*constraint + 1e-4*meandensity*VectorXd::Ones(gd.NxNyNz);
  potential = -potential.cwise().log();
    
  Minimizer min = Precision(1e-6, 
                            PreconditionedConjugateGradient(f, gd, 1, 
                                                            &potential,
                                                            QuadraticLineMinimizer));
    
  double mumax = mu, mumin = mu, dmu = 4.0/N;
  double Nnow = N_from_mu(&min, &potential, constraint, mu);
  const double fraccuracy = 1e-3;
  if (fabs(Nnow/N - 1) > fraccuracy) {
    if (Nnow > N) {
      while (Nnow > N) {
        mumin = mumax;
        mumax += dmu;
        dmu *= 2;
        
        Nnow = N_from_mu(&min, &potential, constraint, mumax);
        // Grid density(gd, EffectivePotentialToDensity()(1, gd, potential));
        // density = EffectivePotentialToDensity()(1, gd, potential);
        // density.epsNativeSlice("papers/contact/figs/box.eps", 
        //                        Cartesian(0,ymax+2,0), Cartesian(0,0,zmax+2), 
        //                        Cartesian(0,-ymax/2-1,-zmax/2-1));
        // density.epsNativeSlice("papers/contact/figs/box-diagonal.eps", 
        //                        Cartesian(xmax+2,0,zmax+2),  Cartesian(0,ymax+2,0),
        //                        Cartesian(-xmax/2-1,-ymax/2-1,-zmax/2-1));
        printf("mumax %g gives N %g\n", mumax, Nnow);
        took("Finding N from mu");
      }
      printf("mu is between %g and %g\n", mumin, mumax);
    } else {
      while (Nnow < N) {
        mumax = mumin;
        if (mumin > dmu) {
          mumin -= dmu;
          dmu *= 2;
        } else if (mumin > 0) {
          mumin = -mumin;
        } else {
          mumin *= 2;
        }
        
        Nnow = N_from_mu(&min, &potential, constraint, mumin);
        // density = EffectivePotentialToDensity()(1, gd, potential);
        // density.epsNativeSlice("papers/contact/figs/box.eps", 
        //                        Cartesian(0,ymax+2,0), Cartesian(0,0,zmax+2), 
        //                        Cartesian(0,-ymax/2-1,-zmax/2-1));
        // density.epsNativeSlice("papers/contact/figs/box-diagonal.eps", 
        //                        Cartesian(xmax+2,0,zmax+2),  Cartesian(0,ymax+2,0),
//.........这里部分代码省略.........

  int FitNullModel(Matrix& mat_Xnull, Matrix& mat_y,
                   const EigenMatrix& kinshipU, const EigenMatrix& kinshipS){
    // sanity check
    if (mat_Xnull.rows != mat_y.rows) return -1;
    if (mat_Xnull.rows != kinshipU.mat.rows()) return -1;
    if (mat_Xnull.rows != kinshipS.mat.rows()) return -1;
    
    // type conversion
    G_to_Eigen(mat_Xnull, &this->ux);
    G_to_Eigen(mat_y, &this->uy);
    this->lambda = kinshipS.mat;
    const Eigen::MatrixXf& U = kinshipU.mat;

    // rotate
    this->ux = U.transpose() * this->ux;
    this->uy = U.transpose() * this->uy;

    // get beta, sigma and delta
    // where delta = sigma2_e / sigma2_g
    double loglik[101];
    int maxIndex = -1;
    double maxLogLik = 0;
    for (int i = 0; i <= 100; ++i ){
      delta = exp(-10 + i * 0.2);
      getBetaSigma2(delta);
      loglik[i] = getLogLikelihood(delta);
#ifdef DEBUG
      fprintf(stderr, "%d\tdelta=%g\tll=%lf\n", i, delta, loglik[i]);
      fprintf(stderr, "beta(0)=%lf\tsigma2=%lf\n",
              beta(0), sigma2);
#endif
      if (std::isnan(loglik[i])) {
        continue;
      }
      if (maxIndex < 0 || loglik[i] > maxLogLik) {
        maxIndex = i;
        maxLogLik = loglik[i];
      }
    }
    if (maxIndex < -1) {
      fprintf(stderr, "Cannot optimize\n");
      return -1;
    }
#if 0
    fprintf(stderr, "maxIndex = %d\tll=%lf\t\tbeta(0)=%lf\tsigma2=%lf\n",
            maxIndex, maxLogLik, beta(0), sigma2);
#endif

    if (maxIndex == 0 || maxIndex == 100) {
      // on the boundary
      // do not try maximize it.
    } else {
      gsl_function F;
      F.function = goalFunction;
      F.params = this;

      Minimizer minimizer;
      double lb = exp(-10 + (maxIndex-1) * 0.2);
      double ub = exp(-10 + (maxIndex+1) * 0.2);
      double start =  exp(-10 + maxIndex * 0.2);
      if (minimizer.minimize(F, start, lb, ub)) {
        fprintf(stderr, "Minimization failed, fall back to initial guess.\n");
        this->delta = start;
      } else {
        this->delta = minimizer.getX();
#ifdef DEBUG       
        fprintf(stderr, "minimization succeed when delta = %g, sigma2 = %g\n", this->delta, this->sigma2);
#endif
      }
    }
    // store some intermediate results
#ifdef DEBUG       
    fprintf(stderr, "delta = sigma2_e/sigma2_g, and sigma2 is sigma2_g\n");
    fprintf(stderr, "maxIndex = %d, delta = %g, Try brent\n", maxIndex, delta);
    fprintf(stderr, "beta[0][0] = %g\t sigma2_g = %g\tsigma2_e = %g\n", beta(0,0), this->sigma2, delta * sigma2);
#endif
    // if (this->test == MetaCov::LRT) {
    // this->nullLikelihood = getLogLikelihood(this->delta);
    // } else if (this->test == MetaCov::SCORE) {
    //   this->uResid = this->uy - this->ux * this->beta;
    // }
    return 0;
  }

void run_with_eta(double eta, const char *name, Functional fhs) {
  // Generates a data file for the pair distribution function, for filling fraction eta
  // and distance of first sphere from wall of z0. Data saved in a table such that the
  // columns are x values and rows are z1 values.
  printf("Now starting run_with_eta with eta = %g name = %s\n",
         eta, name);
  Functional f = OfEffectivePotential(fhs + IdealGas());
  double mu = find_chemical_potential(f, 1, eta/(4*M_PI/3));
  f = OfEffectivePotential(fhs + IdealGas()
                           + ChemicalPotential(mu));
  Lattice lat(Cartesian(width,0,0), Cartesian(0,width,0), Cartesian(0,0,width));
  GridDescription gd(lat, dx);
  Grid potential(gd);
  Grid constraint(gd);
  constraint.Set(notinsphere);
  f = constrain(constraint, f);
  potential = (eta*constraint + 1e-4*eta*VectorXd::Ones(gd.NxNyNz))/(4*M_PI/3);
  potential = -potential.cwise().log();

  const double approx_energy = (fhs + IdealGas() + ChemicalPotential(mu))(1, eta/(4*M_PI/3))*uipow(width,3);
  const double precision = fabs(approx_energy*1e-10);
  //printf("Minimizing to %g absolute precision...\n", precision);
  { // Put mimizer in block so as to free it when we finish minimizing to save memory.
    Minimizer min = Precision(precision,
                              PreconditionedConjugateGradient(f, gd, 1,
                                                              &potential,
                                                              QuadraticLineMinimizer));
    for (int i=0;min.improve_energy(true) && i<100;i++) {
      double peak = peak_memory()/1024.0/1024;
      double current = current_memory()/1024.0/1024;
      printf("Peak memory use is %g M (current is %g M)\n", peak, current);
      fflush(stdout);
    }
    took("Doing the minimization");
  }
  Grid density(gd, EffectivePotentialToDensity()(1, gd, potential));
  Grid gsigma(gd, gSigmaA(1.0)(1, gd, density));
  Grid nA(gd, ShellConvolve(2)(1, density)/(4*M_PI*4));
  Grid n3(gd, StepConvolve(1)(1, density));
  Grid nbar_sokolowski(gd, StepConvolve(1.6)(1, density));
  nbar_sokolowski /= (4.0/3.0*M_PI*ipow(1.6, 3));
  // Create the walls directory if it doesn't exist.
  if (mkdir("papers/pair-correlation/figs/walls", 0777) != 0 && errno != EEXIST) {
    // We failed to create the directory, and it doesn't exist.
    printf("Failed to create papers/pair-correlation/figs/walls: %s",
           strerror(errno));
    exit(1); // fail immediately with error code
  }

  // here you choose the values of z0 to use
  // dx is the resolution at which we compute the density.
  char *plotname = new char[4096];
  for (double z0 = 2.1; z0 < 4.5; z0 += 2.1) {
    // For each z0, we now pick one of our methods for computing the
    // pair distribution function:
    for (int version = 0; version < numplots; version++) {
      sprintf(plotname,
              "papers/pair-correlation/figs/triplet%s-%s-%04.2f-%1.2f.dat",
              name, fun[version], eta, z0);
      FILE *out = fopen(plotname,"w");
      FILE *xfile = fopen("papers/pair-correlation/figs/triplet-x.dat","w");
      FILE *zfile = fopen("papers/pair-correlation/figs/triplet-z.dat","w");
      // the +1 for z0 and z1 are to shift the plot over, so that a sphere touching the wall
      // is at z = 0, to match with the monte carlo data
      const Cartesian r0(0,0,z0);
      for (double x = 0; x < 4; x += dx) {
        for (double z1 = -4; z1 <= 9; z1 += dx) {
          const Cartesian r1(x,0,z1);
          double g2 = pairdists[version](gsigma, density, nA, n3, nbar_sokolowski, r0, r1);
          double n_bulk = (3.0/4.0/M_PI)*eta;
          double g3 = g2*density(r0)*density(r1)/n_bulk/n_bulk;
          fprintf(out, "%g\t", g3);
          fprintf(xfile, "%g\t", x);
          fprintf(zfile, "%g\t", z1);
        }
        fprintf(out, "\n");
        fprintf(xfile, "\n");
        fprintf(zfile, "\n");
      }
      fclose(out);
      fclose(xfile);
      fclose(zfile);
    }
  }
  delete[] plotname;
  took("Dumping the triplet dist plots");
  const double ds = 0.01; // step size to use in path plots, FIXME increase for publication!
  const double delta = .1; //this is the value of radius of the
                           //particle as it moves around the contact
                           //sphere on its path
  char *plotname_path = new char[4096];
  for (int version = 0; version < numplots; version++) {
    sprintf(plotname_path,
            "papers/pair-correlation/figs/triplet%s-path-%s-%04.2f.dat",
            name, fun[version], eta);
    FILE *out_path = fopen(plotname_path, "w");
    if (!out_path) {
      fprintf(stderr, "Unable to create file %s!\n", plotname_path);
      return;
    }
//.........这里部分代码省略.........

int main(int argc, char *argv[]) {
  if (argc > 1) {
    if (sscanf(argv[1], "%lg", &diameter) != 1) {
      printf("Got bad argument: %s\n", argv[1]);
      return 1;
    }
    diameter *= nm;
    using_default_diameter = false;
    printf("Diameter is %g bohr\n", diameter);
  }
  
  const double ptransition  =(3.0*M_PI-4.0)*(diameter/2.0)/2.0;
  const double dmax = ptransition + 0.6*nm;
  double zmax = 2*diameter+dmax+2*nm;
  double ymax = 2*diameter+dmax+2*nm;

  char *datname = new char[1024];
  snprintf(datname, 1024, "papers/water-saft/figs/four-rods-in-water-%04.1fnm.dat", diameter/nm);
  
  FILE *o = fopen(datname, "w");
  delete[] datname;

  Functional f = OfEffectivePotential(WaterSaft(new_water_prop.lengthscale,
                                                new_water_prop.epsilonAB, new_water_prop.kappaAB,
                                                new_water_prop.epsilon_dispersion,
                                                new_water_prop.lambda_dispersion,
                                                new_water_prop.length_scaling, 0));
  double n_1atm = pressure_to_density(f, new_water_prop.kT, atmospheric_pressure,
                                      0.001, 0.01);

  double mu_satp = find_chemical_potential(f, new_water_prop.kT, n_1atm);

  f = OfEffectivePotential(WaterSaft(new_water_prop.lengthscale,
                                     new_water_prop.epsilonAB, new_water_prop.kappaAB,
                                     new_water_prop.epsilon_dispersion,
                                     new_water_prop.lambda_dispersion,
                                     new_water_prop.length_scaling, mu_satp));
  
  const double EperVolume = f(new_water_prop.kT, -new_water_prop.kT*log(n_1atm));
  const double EperCell = EperVolume*(zmax*ymax - 4*0.25*M_PI*diameter*diameter)*width;

  //Functional X = Xassociation(new_water_prop.lengthscale, new_water_prop.epsilonAB, 
  //                            new_water_prop.kappaAB, new_water_prop.epsilon_dispersion,
  //                            new_water_prop.lambda_dispersion,
  //                            new_water_prop.length_scaling);
  
  Functional S = OfEffectivePotential(EntropySaftFluid2(new_water_prop.lengthscale, 
                                                        new_water_prop.epsilonAB,
                                                        new_water_prop.kappaAB,
                                                        new_water_prop.epsilon_dispersion,
                                                        new_water_prop.lambda_dispersion,
                                                        new_water_prop.length_scaling));
  //dmax, dstep already in bohrs (so it doesn't need to be converted from nm)
  double dstep = 0.25*nm;
  for (distance=0.0*nm; distance<=dmax; distance += dstep) {

    if ((distance >= ptransition - 0.5*nm) && (distance <= ptransition + 0.05*nm)) {
      if (distance >= ptransition - 0.25*nm) {
        dstep = 0.03*nm;
      } else {
        dstep = 0.08*nm;
      }
    } else {
      dstep = 0.25*nm;
    }

    Lattice lat(Cartesian(width,0,0), Cartesian(0,ymax,0), Cartesian(0,0,zmax));
    GridDescription gd(lat, 0.2);
    printf("Grid is %d x %d x %d\n", gd.Nx, gd.Ny, gd.Nz);
    Grid potential(gd);
    Grid constraint(gd);
    constraint.Set(notinwall);
    
    f = OfEffectivePotential(WaterSaft(new_water_prop.lengthscale,
                                       new_water_prop.epsilonAB, new_water_prop.kappaAB,
                                       new_water_prop.epsilon_dispersion,
                                       new_water_prop.lambda_dispersion,
                                       new_water_prop.length_scaling, mu_satp));
    f = constrain(constraint, f);
    
    printf("Diameter is %g bohr (%g nm)\n", diameter, diameter/nm);
    printf("Distance between rods = %g bohr (%g nm)\n", distance, distance/nm);


    potential = new_water_prop.liquid_density*constraint
      + 400*new_water_prop.vapor_density*VectorXd::Ones(gd.NxNyNz);
    //potential = new_water_prop.liquid_density*VectorXd::Ones(gd.NxNyNz);
    potential = -new_water_prop.kT*potential.cwise().log();

    const double surface_tension = 5e-5; // crude guess from memory...
    const double surfprecision = 1e-5*(4*M_PI*diameter)*width*surface_tension; // five digits accuracy
    const double bulkprecision = 1e-12*fabs(EperCell); // but there's a limit on our precision for small rods
    const double precision = bulkprecision + surfprecision;
    printf("Precision limit from surface tension is to %g based on %g and %g\n",
           precision, surfprecision, bulkprecision);
    Minimizer min = Precision(precision, PreconditionedConjugateGradient(f, gd, new_water_prop.kT,
                                                                     &potential,
                                                                     QuadraticLineMinimizer));
    const int numiters = 200;
    for (int i=0;i<numiters && min.improve_energy(false);i++) {
//.........这里部分代码省略.........