Commit 4ea3308f authored by Sikandar Y. Mashayak's avatar Sikandar Y. Mashayak

potentialfunctions: moving potential functions from tools to libcsg

parent 8e775fbd
file(GLOB_RECURSE VOTCA_HEADERS *.h)
file(GLOB_RECURSE VOTCA_HEADERS *.h potentialfunctions/*.h)
install(FILES ${VOTCA_HEADERS} DESTINATION include/votca/csg)
/*
*
/*
*
* Author: sikandar
*
* Created on November 8, 2011, 11:52 PM
......@@ -23,7 +23,7 @@ using namespace votca::tools;
class PotentialFunction {
public:
virtual ~PotentialFunction() {}
// read parameters from the input file
virtual void setParam(string filename);
......@@ -39,7 +39,7 @@ public:
// set ith parameter
void setParam(const int i, const double val) { _lam(i) = val; }
// set ith parameter among those to be optimized
virtual void setOptParam(const int i, const double val) {
virtual void setOptParam(const int i, const double val) {
setParam(i,val);
}
// set minimum r value to avoid large values
......@@ -69,9 +69,9 @@ public:
double getMinDist() const { return _min; }
protected:
PotentialFunction(const string& name_,const int nlam_,const double min_,const double max_);
string _name;
ub::vector<double> _lam;
double _cut_off;
......
......@@ -22,20 +22,20 @@ public:
double CalculateDF(const int i, const double r) const;
// calculate second derivative w.r.t. ith parameter
double CalculateD2F(const int i, const int j, const double r) const;
int getOptParamSize() const ;
void setParam(string filename);
// save parameters to the file
void SaveParam(const string& filename);
void SavePotTab(const string& filename, const double step);
void SavePotTab(const string& filename, const double step, const double rmin, const double rcut);
void SavePotTab(const string& filename, const double step, const double rmin, const double rcut);
void setOptParam(const int i, const double val);
double getOptParam(const int i) const;
void extrapolExclParam();
protected:
......@@ -51,10 +51,10 @@ protected:
int _nbreak;
double _dr;
ub::vector<double> _rbreak;
ub::matrix<double> _M;
string _extrapol_type;
string _extrapol_type;
};
......
/*
/*
* Author: sikandar
*
* Created on November 8, 2011, 11:52 PM
......@@ -28,7 +28,7 @@ public:
double CalculateDF(const int i, const double r) const;
// calculate second derivative w.r.t. ith parameter
double CalculateD2F(const int i, const int j, const double r) const;
};
#endif
......@@ -16,9 +16,9 @@ if (WITH_GMX)
endif(GMX_DOUBLE)
if (WITH_GMX_DEVEL)
message(WARNING " Using gromacs-5.0, this is very experimental, hope you know what you are doing")
set(LIBGMX "libgromacs")
set(LIBGMX "libgromacs")
else(WITH_GMX_DEVEL)
set(LIBGMX "libgmx")
set(LIBGMX "libgmx")
endif(WITH_GMX_DEVEL)
find_package(GROMACS COMPONENTS "${LIBGMX}${GMX_SUFFIX}")
if (NOT GROMACS_FOUND)
......@@ -30,7 +30,7 @@ if (WITH_GMX)
set(GMX ${GROMACS_VERSION})
add_executable(gmx_print_version modules/io/gmx_print_version.cc)
target_link_libraries(gmx_print_version ${GROMACS_LIBRARIES})
add_custom_command(OUTPUT gmx_libs_version.h COMMAND gmx_print_version > gmx_libs_version.h DEPENDS gmx_print_version)
add_custom_command(OUTPUT gmx_libs_version.h COMMAND gmx_print_version > gmx_libs_version.h DEPENDS gmx_print_version)
list(APPEND GMX_SOURCES gmx_libs_version.h)
else(WITH_GMX)
set(GMX)
......@@ -43,7 +43,7 @@ include_directories(${CMAKE_CURRENT_BINARY_DIR})
add_custom_target(hgversion COMMAND ${CMAKE_COMMAND} -DTOP_SOURCE_DIR="${CMAKE_SOURCE_DIR}" -P ${CMAKE_MODULE_PATH}/hgversion.cmake)
set_property(DIRECTORY APPEND PROPERTY ADDITIONAL_MAKE_CLEAN_FILES hgversion.h)
file(GLOB VOTCA_SOURCES *.cc)
file(GLOB VOTCA_SOURCES *.cc potentialfunctions/*.cc)
file(GLOB NOT_VOTCA_SOURCES version_nb.cc test.cc)
list(REMOVE_ITEM VOTCA_SOURCES ${NOT_VOTCA_SOURCES})
add_library(votca_csg ${VOTCA_SOURCES} ${GMX_SOURCES} ${IO_SOURCES})
......
#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;
char flag;
for (r_init = rmin, i = 0; i < ngrid - 1; r_init += step) {
// put point, result, flag at point out_x into the table
if( r_init < _min || r_init > _cut_off )
flag = 'o';
else
flag = 'i';
pot_tab.set(i++, r_init, CalculateF(r_init), flag);
}
if( rcut < _min || rcut > _cut_off )
flag = 'o';
else
flag = 'i';
pot_tab.set(i, rcut, CalculateF(rcut), flag);
// save table in the file
pot_tab.Save(filename);
}
#include <votca/csg/potentialfunctions/potentialfunctioncbspl.h>
PotentialFunctionCBSPL::PotentialFunctionCBSPL(const string& name_,const int nlam_,
const int ncutcoeff_, const string& extrapol_type_,
const double min_, const double max_) :
PotentialFunction(name_,nlam_,min_,max_) {
/* Here nlam_ is the total number of coeff values that are to be optimized
* To ensure that potential and force go to zero smoothly near cut-off,
* as suggested in Ref. PCCP, 11, 1901, 2009, coeff values leading up to
* cut-off and beyond take a value of zero.
*
* Since region less than rmin is not sampled sufficiently for stability
* first _nexcl coefficients are not optimized instead their values are
* extrapolated from first statistically significant knot values near rmin
*/
_extrapol_type = extrapol_type_;
// number of break points = _lam.size() - 2
_nbreak = _lam.size() - 2;
_dr = ( _cut_off )/( double (_nbreak - 1) );
// break point locations
// since ncoeff = nbreak +2 , r values for last two coefficients are also
// computed
_rbreak.resize(_lam.size(),false);
_rbreak.clear();
for( int i = 0; i < _lam.size(); i++)
_rbreak(i) = i*_dr;
// exclude knots corresponding to r <= _min
_nexcl = min( int( ( _min )/_dr ), _nbreak - 2 ) + 1;
// account for finite numerical division of _min/_dr
// e.g. 0.24/0.02 may result in 11.99999999999999
if( _rbreak(_nexcl) == _min ) _nexcl++;
_ncutcoeff = ncutcoeff_;
_M.resize(4,4,false);
_M.clear();
_M(0,0) = 1.0; _M(0,1) = 4.0; _M(0,2) = 1.0; _M(0,3) = 0.0;
_M(1,0) = -3.0; _M(1,1) = 0.0; _M(1,2) = 3.0; _M(1,3) = 0.0;
_M(2,0) = 3.0; _M(2,1) = -6.0; _M(2,2) = 3.0; _M(2,3) = 0.0;
_M(3,0) = -1.0; _M(3,1) = 3.0; _M(3,2) = -3.0; _M(3,3) = 1.0;
_M /= 6.0;
}
int PotentialFunctionCBSPL::getOptParamSize() const {
return _lam.size() - _nexcl - _ncutcoeff;
}
void PotentialFunctionCBSPL::setParam(string filename) {
Table param;
param.Load(filename);
_lam.clear();
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 {
// force last _ncutcoeff to be zero
for( int i = 0; i < _lam.size() - _ncutcoeff; i++){
_rbreak(i) = param.x(i);
_lam(i) = param.y(i);
}
}
}
void PotentialFunctionCBSPL::SaveParam(const string& filename){
extrapolExclParam();
Table param;
param.SetHasYErr(false);
param.resize(_lam.size(), false);
// write extrapolated knots with flag 'o'
for (int i = 0; i < _nexcl; i++)
param.set(i, _rbreak(i), _lam(i), 'o');
for (int i = _nexcl; i < _lam.size(); i++)
param.set(i, _rbreak(i), _lam(i), 'i');
param.Save(filename);
}
void PotentialFunctionCBSPL::SavePotTab(const string& filename,
const double step, const double rmin, const double rcut) {
extrapolExclParam();
PotentialFunction::SavePotTab(filename,step,rmin,rcut);
}
void PotentialFunctionCBSPL::SavePotTab(const string& filename,
const double step) {
extrapolExclParam();
PotentialFunction::SavePotTab(filename,step);
}
void PotentialFunctionCBSPL::extrapolExclParam(){
if( _extrapol_type == "linear") {
// extrapolate first _nexcl knot values using linear extrapolation
// u(r) = ar + b
// a = m
// b = - m*r0 + u0
// m = (u1-u0)/(r1-r0)
double u0 = _lam(_nexcl);
double r0 = _rbreak(_nexcl);
double m = (_lam(_nexcl + 1) - _lam(_nexcl)) /
(_rbreak(_nexcl + 1) - _rbreak(_nexcl));
// check for positive slope
if( m >= 0.0 )
throw std::runtime_error("In potential "+_name+": min r value for cbspl is too large,\n"
"for linear extrapolation, choose min r such that potential slope at min < 0,\n"
"else extrapolated knot values in "
"the repulsive core would be negative or zero.\n");
double a = m;
double b = -1.0*m*r0 + u0;
for (int i = 0; i < _nexcl; i++)
_lam(i) = a*_rbreak(i) + b;
} else if (_extrapol_type=="exponential"){
// extrapolate first _nexcl knot values using exponential extrapolation
// u(r) = a * exp( b * r)
// a = u0 * exp ( - m * r0/u0 )
// b = m/u0
// m = (u1-u0)/(r1-r0)
double u0 = _lam(_nexcl);
if( u0 <= 0.0 )
throw std::runtime_error("In potential "+_name+": min r value for cbspl is too large,\n"
"for exponential extrapolation, choose min r such that knot value at (rmin+dr) > 0,\n"
"else extrapolated knot values in "
"the repulsive core would be negative or zero.\n");
double r0 = _rbreak(_nexcl);
double m = (_lam(_nexcl + 1) - _lam(_nexcl)) /
(_rbreak(_nexcl + 1) - _rbreak(_nexcl));
double a = u0 * exp(-m * r0 / u0);
double b = m / u0;
for (int i = 0; i < _nexcl; i++)
_lam(i) = a * exp(b * _rbreak(i));
} else
throw std::runtime_error("Extrapolation method \""
+ _extrapol_type + "\" selected for \""
+ _name + "\" is not available yet.\n"
+ "Please specify either \"linear or exponential\" "
+ " in options file.");
}
void PotentialFunctionCBSPL::setOptParam(const int i, const double val){
_lam( i + _nexcl ) = val;
}
double PotentialFunctionCBSPL::getOptParam(const int i) const{
return _lam( i + _nexcl );
}
double PotentialFunctionCBSPL::CalculateF (const double r) const {
if( r <= _cut_off){
ub::vector<double> R;
ub::vector<double> B;
R.resize(4,false); R.clear();
B.resize(4,false); B.clear();
int indx = min( int( r /_dr ) , _nbreak-2 );
double rk = indx*_dr;
double t = ( r - rk)/_dr;
R(0) = 1.0; R(1) = t; R(2) = t*t; R(3) = t*t*t;
ub::vector<double> RM = ub::prod(R,_M);
B(0) = _lam(indx); B(1) = _lam(indx+1); B(2) = _lam(indx+2);
B(3) = _lam(indx+3);
return ub::inner_prod(B,RM);
} else
return 0.0;
}
// calculate first derivative w.r.t. ith parameter
double PotentialFunctionCBSPL::CalculateDF(const int i, const double r) const{
// since first _nexcl parameters are not optimized for stability reasons
//i = i + _nexcl;
if ( r <= _cut_off ) {
int indx = min( int( ( r )/_dr ), _nbreak-2 );
if ( i + _nexcl >= indx && i + _nexcl <= indx+3 ){
ub::vector<double> R;
R.resize(4,false); R.clear();
double rk = indx*_dr;
double t = ( r - rk)/_dr;
R(0) = 1.0; R(1) = t; R(2) = t*t; R(3) = t*t*t;
ub::vector<double> RM = ub::prod(R,_M);
return RM(i + _nexcl-indx);
}else
return 0.0;
} else
return 0;
}
// calculate second derivative w.r.t. ith parameter
double PotentialFunctionCBSPL::CalculateD2F(const int i, const int j,
const double r) const {
// for cubic B-SPlines D2F is zero for all lamdas
return 0.0;
}
#include "potentialfunctionlj126.h"
#include <votca/csg/potentialfunctions/potentialfunctionlj126.h>
PotentialFunctionLJ126::PotentialFunctionLJ126(const string& name_,const double min_,
const double max_) : PotentialFunction(name_,2,min_,max_){
const double max_) : PotentialFunction(name_,2,min_,max_){
}
double PotentialFunctionLJ126::CalculateF (const double r) const {
if ( r >= _min && r <= _cut_off ) {
return _lam(0)/pow(r,12) - _lam(1)/pow(r,6) ;
} else {
return 0.0;
if ( r >= _min && r <= _cut_off )
return _lam(0)/pow(r,12) - _lam(1)/pow(r,6) ;
else
return 0.0;
}
}
// calculate first derivative w.r.t. ith parameter
double PotentialFunctionLJ126::CalculateDF(const int i, const double r) const {
if ( r >= _min && r <= _cut_off ) {
if ( r >= _min && r <= _cut_off ) {
switch(i) {
case 0:
return 1.0/pow(r,12);
case 1:
return -1.0/pow(r,6);
}
} else {
switch(i) {
case 0:
return 1.0/pow(r,12);
case 1:
return -1.0/pow(r,6);
}
return 0.0;
} else
return 0.0;
}
}
// calculate second derivative w.r.t. ith and jth parameters
double PotentialFunctionLJ126::CalculateD2F(const int i, const int j,
const double r) const {
double PotentialFunctionLJ126::CalculateD2F(const int i, const int j,
const double r) const {
return 0.0;
}
#include <votca/csg/potentialfunctions/potentialfunctionljg.h>
PotentialFunctionLJG::PotentialFunctionLJG(const string& name_,const double min_,
const double max_) : PotentialFunction(name_,5,min_,max_){
}
double PotentialFunctionLJG::CalculateF (const double r) const {
// lj 12-6 part + gaussian
if ( r >= _min && r <= _cut_off )
return _lam(0)/pow(r,12) - _lam(1)/pow(r,6)
+ _lam(2)*exp( -1.0*_lam(3)*(r-_lam(4))*(r-_lam(4)) );
else
return 0.0;
}
// calculate first derivative w.r.t. ith parameter
double PotentialFunctionLJG::CalculateDF(const int i, const double r) const {
if ( r >= _min && r <= _cut_off ) {
switch(i) {
case 0:
return 1.0/pow(r,12);
case 1:
return -1.0/pow(r,6);
case 2:
return exp( -1.0*_lam(3)*(r-_lam(4))*(r-_lam(4)) );
case 3:
return -1.0*_lam(2)*(r-_lam(4))*(r-_lam(4)) *
exp( -1.0*_lam(3)*(r-_lam(4))*(r-_lam(4)) );
case 4:
return 2.0*_lam(2)*_lam(3)*(r-_lam(4)) *
exp( -1.0*_lam(3)*(r-_lam(4))*(r-_lam(4)) );
}
} else
return 0.0;
}
// calculate second derivative w.r.t. ith parameter
double PotentialFunctionLJG::CalculateD2F(const int i, const int j,
const double r) const {
if ( r >= _min && r <= _cut_off ) {
switch(i) {
case 0:
// all second derivatives w.r.t. c12 are zero
return 0.0;
case 1:
return 0.0;
case 2:
switch(j){
case 0:
return 0.0;
case 1:
return 0.0;
case 2:
return 0.0;
case 3:
return -1.0*(r-_lam(4))*(r-_lam(4)) *
exp( -1.0*_lam(3)*(r-_lam(4))*(r-_lam(4)) );
case 4:
return 2.0*_lam(3)*(r-_lam(4)) *
exp( -1.0*_lam(3)*(r-_lam(4))*(r-_lam(4)) );
}
case 3:
switch(j){
case 0:
return 0.0;
case 1:
return 0.0;
case 2:
return -1.0*(r-_lam(4))*(r-_lam(4)) *
exp( -1.0*_lam(3)*(r-_lam(4))*(r-_lam(4)) );
case 3:
return _lam(2)*pow((r-_lam(4)),4) *
exp( -1.0*_lam(3)*(r-_lam(4))*(r-_lam(4)) );
case 4:
return 2.0*_lam(2)*(r-_lam(4)) *
( 1.0 - _lam(3)*pow((r-_lam(4)),2) ) *
exp( -1.0*_lam(3)*(r-_lam(4))*(r-_lam(4)) );
}
case 4:
switch(j){
case 0:
return 0.0;
case 1:
return 0.0;
case 2:
return 2.0*_lam(3)*(r-_lam(4)) *
exp( -1.0*_lam(3)*(r-_lam(4))*(r-_lam(4)) );
case 3:
return 2.0*_lam(2)*(r-_lam(4)) *
( 1.0 - _lam(3)*pow((r-_lam(4)),2) ) *
exp( -1.0*_lam(3)*(r-_lam(4))*(r-_lam(4)) );
case 4:
return 2.0*_lam(2)*_lam(3)*
( 2.0*_lam(3)*pow((r-_lam(4)),2) - 1.0 ) *
exp( -1.0*_lam(3)*(r-_lam(4))*(r-_lam(4)) );
}
}
} else
return 0.0;
}
foreach(PROG csg_reupdate csg_map csg_dump csg_property csg_resample csg_stat csg_fmatch csg_imcrepack csg_gmxtopol csg_density)
file(GLOB ${PROG}_SOURCES ${PROG}*.cc)
file(GLOB csg_reupdate_SOURCES csg_reupdate.cc potentialfunctions/*.cc)
add_executable(${PROG} ${${PROG}_SOURCES})
target_link_libraries(${PROG} votca_csg)
install(TARGETS ${PROG} RUNTIME DESTINATION bin)
if (TXT2TAGS_FOUND AND BASH)
add_custom_command(OUTPUT ${PROG}.man
add_custom_command(OUTPUT ${PROG}.man
COMMAND ${CMAKE_CURRENT_BINARY_DIR}/${PROG} --help > ${PROG}.help
COMMAND ${BASH} ${CMAKE_BINARY_DIR}/scripts/help2t2t ${PROG}.help > ${PROG}.t2t
COMMAND ${BASH} ${CMAKE_BINARY_DIR}/scripts/help2t2t ${PROG}.help > ${PROG}.t2t
COMMAND ${TXT2TAGS_EXECUTABLE} -q -t man -i ${PROG}.t2t -o ${PROG}.man
DEPENDS scripts/help2t2t ${PROG})
add_custom_target(${PROG}_manpage DEPENDS ${PROG}.man)
......