Skip to content

Commit

Permalink
applied new style to tools
Browse files Browse the repository at this point in the history
  • Loading branch information
@Iximiel
Iximiel committed Sep 30, 2024
1 parent cc60f37 commit e9e421b
Show file tree
Hide file tree
Showing 55 changed files with 3,600 additions and 1,513 deletions.
1 change: 0 additions & 1 deletion src/astyle.sh
View file
Original file line number Diff line number Diff line change
Expand Up @@ -40,7 +40,6 @@ use_legacy=(
s2cm
secondarystructure
setup
tools
vatom
ves
vesselbase
Expand Down
97 changes: 70 additions & 27 deletions src/tools/BiasRepresentation.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -92,14 +92,23 @@ BiasRepresentation::BiasRepresentation(const std::vector<Value*> & tmpvalues, Co
void BiasRepresentation::addGrid(const std::vector<std::string> & gmin, const std::vector<std::string> & gmax, const std::vector<unsigned> & nbin ) {
plumed_massert(hills.size()==0,"you can set the grid before loading the hills");
plumed_massert(hasgrid==false,"to build the grid you should not having the grid in this bias representation");
std::string ss; ss="file.free";
std::vector<Value*> vv; for(unsigned i=0; i<values.size(); i++) vv.push_back(values[i]);
std::string ss;
ss="file.free";
std::vector<Value*> vv;
for(unsigned i=0; i<values.size(); i++)
vv.push_back(values[i]);
BiasGrid_=Tools::make_unique<Grid>(ss,vv,gmin,gmax,nbin,false,true);
hasgrid=true;
}

bool BiasRepresentation::hasSigmaInInput() {
if(histosigma.size()==0) {return false;} else {return true;}
//return histosigma.size()!=0;?
if(histosigma.size()==0) {
return false;
}
else {
return true;
}
}

void BiasRepresentation::setRescaledToBias(bool rescaled) {
Expand Down Expand Up @@ -153,14 +162,19 @@ void BiasRepresentation::pushKernel( IFile *ifile ) {
}
// the bias factor is not something about the kernels but
// must be stored to keep the bias/free energy duality
std::string dummy; double dummyd;
std::string dummy;
double dummyd;
if(ifile->FieldExist("biasf")) {
ifile->scanField("biasf",dummy);
Tools::convert(dummy,dummyd);
} else {dummyd=1.0;}
} else {
dummyd=1.0;
}
biasf.push_back(dummyd);
// the domain does not pertain to the kernel but to the values here defined
std::string mins,maxs,minv,maxv,mini,maxi; mins="min_"; maxs="max_";
std::string mins,maxs,minv,maxv,mini,maxi;
mins="min_";
maxs="max_";
for(int i=0 ; i<ndim; i++) {
if(values[i]->isPeriodic()) {
ifile->scanField(mins+names[i],minv);
Expand All @@ -177,24 +191,33 @@ void BiasRepresentation::pushKernel( IFile *ifile ) {
std::vector<unsigned> nneighb;
if(doInt_&&(kk->getCenter()[0]+kk->getContinuousSupport()[0] > uppI_ || kk->getCenter()[0]-kk->getContinuousSupport()[0] < lowI_ )) {
nneighb=BiasGrid_->getNbin();
} else nneighb=kk->getSupport(BiasGrid_->getDx());
} else
nneighb=kk->getSupport(BiasGrid_->getDx());
std::vector<Grid::index_t> neighbors=BiasGrid_->getNeighbors(kk->getCenter(),nneighb);
std::vector<double> der(ndim);
std::vector<double> xx(ndim);
if(mycomm.Get_size()==1) {
for(unsigned i=0; i<neighbors.size(); ++i) {
Grid::index_t ineigh=neighbors[i];
for(int j=0; j<ndim; ++j) {der[j]=0.0;}
for(int j=0; j<ndim; ++j) {
der[j]=0.0;
}
BiasGrid_->getPoint(ineigh,xx);
// assign xx to a new vector of values
for(int j=0; j<ndim; ++j) {values[j]->set(xx[j]);}
for(int j=0; j<ndim; ++j) {
values[j]->set(xx[j]);
}
double bias;
if(doInt_) bias=kk->evaluate(values,der,true,doInt_,lowI_,uppI_);
else bias=kk->evaluate(values,der,true);
if(doInt_)
bias=kk->evaluate(values,der,true,doInt_,lowI_,uppI_);
else
bias=kk->evaluate(values,der,true);
if(rescaledToBias) {
double f=(biasf.back()-1.)/(biasf.back());
bias*=f;
for(int j=0; j<ndim; ++j) {der[j]*=f;}
for(int j=0; j<ndim; ++j) {
der[j]*=f;
}
}
BiasGrid_->addValueAndDerivatives(ineigh,bias,der);
}
Expand All @@ -207,23 +230,34 @@ void BiasRepresentation::pushKernel( IFile *ifile ) {
for(unsigned i=rank; i<neighbors.size(); i+=stride) {
Grid::index_t ineigh=neighbors[i];
BiasGrid_->getPoint(ineigh,xx);
for(int j=0; j<ndim; ++j) {values[j]->set(xx[j]);}
if(doInt_) allbias[i]=kk->evaluate(values,der,true,doInt_,lowI_,uppI_);
else allbias[i]=kk->evaluate(values,der,true);
for(int j=0; j<ndim; ++j) {
values[j]->set(xx[j]);
}
if(doInt_)
allbias[i]=kk->evaluate(values,der,true,doInt_,lowI_,uppI_);
else
allbias[i]=kk->evaluate(values,der,true);
if(rescaledToBias) {
double f=(biasf.back()-1.)/(biasf.back());
allbias[i]*=f;
for(int j=0; j<ndim; ++j) {tmpder[j]*=f;}
for(int j=0; j<ndim; ++j) {
tmpder[j]*=f;
}
}
// this solution with the temporary vector is rather bad, probably better to take
// a pointer of double as it was in old gaussian
for(int j=0; j<ndim; ++j) { allder[ndim*i+j]=tmpder[j]; tmpder[j]=0.;}
for(int j=0; j<ndim; ++j) {
allder[ndim*i+j]=tmpder[j];
tmpder[j]=0.;
}
}
mycomm.Sum(allbias);
mycomm.Sum(allder);
for(unsigned i=0; i<neighbors.size(); ++i) {
Grid::index_t ineigh=neighbors[i];
for(int j=0; j<ndim; ++j) {der[j]=allder[ndim*i+j];}
for(int j=0; j<ndim; ++j) {
der[j]=allder[ndim*i+j];
}
BiasGrid_->addValueAndDerivatives(ineigh,allbias[i],der);
}
}
Expand All @@ -242,10 +276,14 @@ Grid* BiasRepresentation::getGridPtr() {

void BiasRepresentation::getMinMaxBin(std::vector<double> &vmin, std::vector<double> &vmax, std::vector<unsigned> &vbin) {
std::vector<double> ss,cc,binsize;
vmin.clear(); vmin.resize(ndim,10.e20);
vmax.clear(); vmax.resize(ndim,-10.e20);
vbin.clear(); vbin.resize(ndim);
binsize.clear(); binsize.resize(ndim,10.e20);
vmin.clear();
vmin.resize(ndim,10.e20);
vmax.clear();
vmax.resize(ndim,-10.e20);
vbin.clear();
vbin.resize(ndim);
binsize.clear();
binsize.resize(ndim,10.e20);
int ndiv=10; // adjustable parameter: division per support
for(unsigned i=0; i<hills.size(); i++) {
if(histosigma.size()!=0) {
Expand All @@ -258,18 +296,23 @@ void BiasRepresentation::getMinMaxBin(std::vector<double> &vmin, std::vector<dou
double dmin=cc[j]-ss[j];
double dmax=cc[j]+ss[j];
double ddiv=ss[j]/double(ndiv);
if(dmin<vmin[j])vmin[j]=dmin;
if(dmax>vmax[j])vmax[j]=dmax;
if(ddiv<binsize[j])binsize[j]=ddiv;
if(dmin<vmin[j])
vmin[j]=dmin;
if(dmax>vmax[j])
vmax[j]=dmax;
if(ddiv<binsize[j])
binsize[j]=ddiv;
}
}
for(int j=0; j<ndim; j++) {
// reset to periodicity
if(values[j]->isPeriodic()) {
double minv,maxv;
values[j]->getDomain(minv,maxv);
if(minv>vmin[j])vmin[j]=minv;
if(maxv<vmax[j])vmax[j]=maxv;
if(minv>vmin[j])
vmin[j]=minv;
if(maxv<vmax[j])
vmax[j]=maxv;
}
vbin[j]=static_cast<unsigned>(std::ceil((vmax[j]-vmin[j])/binsize[j]) );
}
Expand Down
66 changes: 49 additions & 17 deletions src/tools/Brent1DRootSearch.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -74,8 +74,13 @@ Brent1DRootSearch<FCLASS>::Brent1DRootSearch( const FCLASS& pf, const double& t

template <class FCLASS>
void Brent1DRootSearch<FCLASS>::bracket( const double& a, const double& b, eng_pointer eng ) {
plumed_assert( a!=b ); ax=a; bx=b; fa=(myclass_func.*eng)(a); fb=(myclass_func.*eng)(b);
if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0)) plumed_merror("input points do not bracket root");
plumed_assert( a!=b );
ax=a;
bx=b;
fa=(myclass_func.*eng)(a);
fb=(myclass_func.*eng)(b);
if ((fa > 0.0 && fb > 0.0) || (fa < 0.0 && fb < 0.0))
plumed_merror("input points do not bracket root");
bracketed=true;
}

Expand All @@ -85,31 +90,58 @@ double Brent1DRootSearch<FCLASS>::search( eng_pointer eng ) {

double cx=bx, d, e, min1, min2, fc=fb, p, q, r, s, tol1, xm;
for(unsigned iter=0; iter<ITMAX; iter++) {
if ( (fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0) ) { cx=ax; fc=fa; e=d=bx-ax; }
if( std::fabs(fc) < std::fabs(fb) ) { ax=bx; bx=cx; cx=ax; fa=fb; fb=fc; fc=fa; }
tol1=2*EPS*std::fabs(bx)+0.5*tol; xm=0.5*(cx-bx);
if( std::fabs(xm) <= tol1 || fb == 0.0 ) return bx;
if ( (fb>0.0 && fc>0.0) || (fb<0.0 && fc<0.0) ) {
cx=ax;
fc=fa;
e=d=bx-ax;
}
if( std::fabs(fc) < std::fabs(fb) ) {
ax=bx;
bx=cx;
cx=ax;
fa=fb;
fb=fc;
fc=fa;
}
tol1=2*EPS*std::fabs(bx)+0.5*tol;
xm=0.5*(cx-bx);
if( std::fabs(xm) <= tol1 || fb == 0.0 )
return bx;
if( std::fabs(e) >= tol1 && std::fabs(fa) > std::fabs(fb) ) {
s=fb/fa;
if( ax==cx ) {
p=2.0*xm*s; q=1.0-s;
p=2.0*xm*s;
q=1.0-s;
} else {
q=fa/fc; r=fb/fc; p=s*(2.0*xm*q*(q-r)-(bx-ax)*(r-1.0)); q=(q-1.0)*(r-1.0)*(s-1.0);
q=fa/fc;
r=fb/fc;
p=s*(2.0*xm*q*(q-r)-(bx-ax)*(r-1.0));
q=(q-1.0)*(r-1.0)*(s-1.0);
}
if (p > 0.0) q = -q;
p=std::fabs(p); min1=3.0*xm*q-std::fabs(tol1*q); min2=std::fabs(e*q);
if (p > 0.0)
q = -q;
p=std::fabs(p);
min1=3.0*xm*q-std::fabs(tol1*q);
min2=std::fabs(e*q);
if (2.0*p < (min1 < min2 ? min1 : min2)) {
e=d; d=p/q;
e=d;
d=p/q;
} else {
d=xm; e=d;
d=xm;
e=d;
}
} else {
d=xm; e=d;
d=xm;
e=d;
}
ax=bx; fa=fb;
if( std::fabs(d) > tol1 ) bx+=d;
else if(xm<0 ) bx -= std::fabs(tol1); // SIGN(tol1,xm);
else bx += tol1;
ax=bx;
fa=fb;
if( std::fabs(d) > tol1 )
bx+=d;
else if(xm<0 )
bx -= std::fabs(tol1); // SIGN(tol1,xm);
else
bx += tol1;
fb = (myclass_func.*eng)(bx);
}

Expand Down
7 changes: 5 additions & 2 deletions src/tools/Citations.cpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -28,8 +28,11 @@ namespace PLMD {

std::string Citations::cite(const std::string & item) {
unsigned i;
for(i=0; i<items.size(); ++i) if(items[i]==item) break;
if(i==items.size()) items.push_back(item);
for(i=0; i<items.size(); ++i)
if(items[i]==item)
break;
if(i==items.size())
items.push_back(item);
plumed_assert(i<items.size());
std::string ret;
Tools::convert(i+1,ret);
Expand Down
Loading

0 comments on commit e9e421b

Please sign in to comment.