Setting Potentials

Here it is explained how to add or modify a potential. The procedure is the same both in the case of a one-particle potential or a interaction term. Lets suppose we want to add a time-independent potential term called mexicanhat. Lets go through the steps needed to implement it:

1. Add the function prototype to the header potential.h. It must take as arguments the parameters entering the potential as type double* and the position of the particle as remaining double variable (in the case of a two-particles interaction potential, there will be just two of them).

   // this is a part of potential.h
   double V1nopotential(double*,double);
   double V1harmonic(double*,double);
   double V1mexicanhat(double*,double); // here the new potential

2. Define the function in the library potential.cc, just below the other. Notice that in this version there is a maximum number of parameters that can be set in the PIMCParams structure (now it is set to 9).

   // this is a part of potential.cc
   double V1nopotential(double*, double){
      return 0.;
   }
   
   double V1harmonic(double* par, double x){
      return par[0]*(x-par[1])*(x-par[1]);
   }
   
   double V1mexicanhat(double* par, double x){  // here the new potential
      return par[0]*(x**2-par[1])**2
   }

3. Add the potential to the wrapper function in potentials.cc, so that the program will bind it and use it in the PIMCClass if its name is indicated in the init file. Notice that there are two wrapper function corresponding, respectively, to the bind of the one-particle potential or of the interaction term.

   // this is a part of ``potential.cc`` corresponding to the V1 wrapper
   std::function<double(double)> V1wrapper(std::string V1form, double* V1par){
   
      using namespace std::placeholders; //to use _1,_2,...
      std::function<double(double)> V1;
      if(V1form=="nopotential") V1 = std::bind(V1nopotential,V1par,_1);
      else if(V1form=="harmonic") V1 = std::bind(V1harmonic,V1par,_1);
      // here the new potential
      else if(V1form=="mexicanhat") V1 = std::bind(V1mexicanhat,V1par,_1);
      else return NULL;
   
      return V1;
   }

Done! You can call the new potential using the name mexicanhat in the init file, specifying the number of parameters (that are 2 in our case) and their value.