potentialfunction.cc 2.44 KB
Newer Older
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81

#include <votca/csg/potentialfunctions/potentialfunction.h>

PotentialFunction::PotentialFunction(const string& name_,const int nlam_,
                                     const double min_,const double max_){

  _name = name_;
  _lam.resize(nlam_);
  _lam.clear();
  _min = min_;
  _cut_off = max_;

}


void PotentialFunction::setParam(string filename) {

  Table param;
  param.Load(filename);

  if( param.size() != _lam.size()) {

    throw std::runtime_error("In potential " + _name +": parameters size mismatch!\n"
                             "Check input parameter file \""
                             + filename + "\" \nThere should be "
                             + boost::lexical_cast<string>( _lam.size() ) + " parameters");
  } else {
    for( int i = 0; i < _lam.size(); i++)
      _lam(i) = param.y(i);
  }

}

void PotentialFunction::SaveParam(const string& filename){

  Table param;
  param.SetHasYErr(false);
  param.resize(_lam.size(), false);

  for (int i = 0; i < _lam.size(); i++)
    param.set(i, i, _lam(i), 'i');

  param.Save(filename);

}

void PotentialFunction::SavePotTab(const string& filename,
                                   const double step){

  int ngrid = (int) ((_cut_off - _min) / step + 1.00000001);

  Table pot_tab;
  pot_tab.SetHasYErr(false);
  pot_tab.resize(ngrid, false);

  double r_init;
  int i;

  // put point, result, flag at point out_x into the table
  for (r_init = _min, i = 0; i < ngrid - 1; r_init += step)
    pot_tab.set(i++, r_init, CalculateF(r_init), 'i');

  pot_tab.set(i, _cut_off, CalculateF(_cut_off), 'i');

  // save table in the file
  pot_tab.Save(filename);

}

void PotentialFunction::SavePotTab(const string& filename,
                                   const double step,
                                   const double rmin, const double rcut){

  int ngrid = (int) ((rcut - rmin) / step + 1.00000001);

  Table pot_tab;
  pot_tab.SetHasYErr(false);
  pot_tab.resize(ngrid, false);

  double r_init;
  int i;
82
  char flag = 'i';
83 84 85 86 87

  for (r_init = rmin, i = 0; i < ngrid - 1; r_init += step) {

    // put point, result, flag at point out_x into the table

88 89 90 91
    //if( r_init < _min || r_init > _cut_off )
    //flag = 'o';
    //else
    //flag = 'i';
92 93 94 95 96

    pot_tab.set(i++, r_init, CalculateF(r_init), flag);

  }

97 98 99 100
  //  if( rcut < _min || rcut > _cut_off )
  //    flag = 'o';
  //  else
  //    flag = 'i';
101 102 103 104 105 106 107

  pot_tab.set(i, rcut, CalculateF(rcut), flag);

  // save table in the file
  pot_tab.Save(filename);

}