Dear All,
I want to use mcintegrate.h to calculate the effective volume of a molecule.

Here is the function to perform this calculation:

#include "C:\Program Files\Numerical Recipes\NR_C302\code\nr3.h"
#include "C:\Program Files\Numerical Recipes\NR_C302\code/ran.h"
#include "C:\Program Files\Numerical Recipes\NR_C302\code/mcintegrate.h"
double dmdcfgVeffMcInt(int *iErrCode, int nSitesPerMolec, double *sigma,
double *rxPos,double *ryPos,double *rzPos,
double *dxMax,double *dyMax,double *dzMax){
//Purpose: compute the effective volume in nm^3 of a molecule defined by sigma and rPositions.

//Input: sigma and rPositions in nm.

/* Written by Amir Vahid 9/18/2010
*/
// transform Neil's discrete mesh into Monte Carlo integration and check speed/accuracy improvement.
int numSites=nSitesPerMolec;
double CH[111][4];//this makes 111 the maximum number of sites per molecule
//todo: make CH allocation dynamic by using malloc or new
*iErrCode=0;

for(int iSite=0;iSite<numSites;iSite++)
{
//PutMsg("Site %d: %lf %lf %lf %lf\n",iSite,sigma[iSite],rxPos[iSite],ryPos[iSite],rzPos[iSite]);
CH[iSite][0]=sigma[iSite];
CH[iSite][1]=rxPos[iSite];
CH[iSite][2]=ryPos[iSite];
CH[iSite][3]=rzPos[iSite];
}


double maxCoordinates[2][3]={{0,0,0},{0,0,0}};
// Create Correct spanning volume to Envelop Complete Molecule
for (int atom=0;atom<numSites;atom++)
{
for (int count=0;count<3;count++)
{
if ((CH[atom][count+1]-(CH[atom][0]/2))<maxCoordinates[0][count]) //jre should subtract radius, not diameter.
maxCoordinates[0][count]=CH[atom][count+1]-(CH[atom][0]/2);
//jre CH[atom][0] is the diameter of this site
if ((CH[atom][count+1]+(CH[atom][0]/2))>maxCoordinates[1][count])
maxCoordinates[1][count]=CH[atom][count+1]+(CH[atom][0]/2);
}
}
*dxMax=fabs(maxCoordinates[1][0]-maxCoordinates[0][0]);
*dyMax=fabs(maxCoordinates[1][1]-maxCoordinates[0][1]);
*dzMax=fabs(maxCoordinates[1][2]-maxCoordinates[0][2]);
double xMax=fabs(maxCoordinates[1][0]-maxCoordinates[0][0]);
double yMax=fabs(maxCoordinates[1][1]-maxCoordinates[0][1]);
double zMax=fabs(maxCoordinates[1][2]-maxCoordinates[0][2]);
VecDoub xlo(3),xhi(3);
xlo[0]=maxCoordinates[0][0];xhi[0]=maxCoordinates[1][0];
xlo[1]=maxCoordinates[0][1];xhi[1]=maxCoordinates[1][1];
xlo[2]=maxCoordinates[0][2];xhi[2]=maxCoordinates[1][2];
const VecDoub x(3);
MCintegrate mymc(xlo,xhi,moleculesfunc,moleculesRegion(x,numSi tes,CH),NULL,10201);
mymc.step(1000000);
mymc.calcanswers();
double vEffNm3=mymc.ff[0];
return vEffNm3;
}

I also have modified the mcintegrate.h as follows:

struct MCintegrate {
Int ndim,nfun,n;
VecDoub ff,fferr;
VecDoub xlo,xhi,x,xx,fn,sf,sferr;
Doub vol;

VecDoub (*funcsp)(const VecDoub &);
VecDoub (*xmapp)(const VecDoub &);
Bool (*inregionp)(const VecDoub &);
Bool (*inregionp_vEff)(const VecDoub &,Int,const MatDoub &);
Ran ran;

MCintegrate(const VecDoub &xlow, const VecDoub &xhigh,
VecDoub funcs(const VecDoub &), Bool inregion(const VecDoub &),
VecDoub xmap(const VecDoub &), Int ranseed);

MCintegrate(const VecDoub &xlow, const VecDoub &xhigh, //overloaded constructor
VecDoub funcs(const VecDoub &), Bool inregion(const VecDoub &,Int,const MatDoub &),
VecDoub xmap(const VecDoub &), Int ranseed);

void step(Int nstep);

void calcanswers();
};
MCintegrate::MCintegrate(const VecDoub &xlow, const VecDoub &xhigh,
VecDoub funcs(const VecDoub &), Bool inregion(const VecDoub &),
VecDoub xmap(const VecDoub &), Int ranseed)
: ndim(xlow.size()), n(0), xlo(xlow), xhi(xhigh), x(ndim), xx(ndim),
funcsp(funcs), xmapp(xmap), inregionp(inregion), vol(1.), ran(ranseed) {
if (xmapp) nfun = funcs(xmapp(xlo)).size();
else nfun = funcs(xlo).size();
ff.resize(nfun);
fferr.resize(nfun);
fn.resize(nfun);
sf.assign(nfun,0.);
sferr.assign(nfun,0.);
for (Int j=0;j<ndim;j++) vol *= abs(xhi[j]-xlo[j]);
}

//AV overloaded constructor
MCintegrate::MCintegrate(const VecDoub &xlow, const VecDoub &xhigh,
VecDoub funcs(const VecDoub &), Bool inregion(const VecDoub &,int,const MatDoub &),
VecDoub xmap(const VecDoub &), Int ranseed)
: ndim(xlow.size()), n(0), xlo(xlow), xhi(xhigh), x(ndim), xx(ndim),
funcsp(funcs), xmapp(xmap), inregionp_vEff(inregion), vol(1.), ran(ranseed) {
if (xmapp) nfun = funcs(xmapp(xlo)).size();
else nfun = funcs(xlo).size();
ff.resize(nfun);
fferr.resize(nfun);
fn.resize(nfun);
sf.assign(nfun,0.);
sferr.assign(nfun,0.);
for (Int j=0;j<ndim;j++) vol *= abs(xhi[j]-xlo[j]);
}


void MCintegrate::step(Int nstep) {
Int i,j;
for (i=0;i<nstep;i++) {
for (j=0;j<ndim;j++)
x[j] = xlo[j]+(xhi[j]-xlo[j])*ran.doub();
if (xmapp) xx = (*xmapp)(x);
else xx = x;
if ((*inregionp)(xx)) {
fn = (*funcsp)(xx);
for (j=0;j<nfun;j++) {
sf[j] += fn[j];
sferr[j] += SQR(fn[j]);
}
}
}
n += nstep;
}

void MCintegrate::calcanswers(){
for (Int j=0;j<nfun;j++) {
ff[j] = vol*sf[j]/n;
fferr[j] = vol*sqrt((sferr[j]/n-SQR(sf[j]/n))/n);
}
}
VecDoub torusfuncs(const VecDoub &x) {
Doub den = 1.;
VecDoub f(4);
f[0] = den;
for (Int i=1;i<4;i++) f[i] = x[i-1]*den;
return f;
}

Bool torusregion(const VecDoub &x) {
return SQR(x[2])+SQR(sqrt(SQR(x[0])+SQR(x[1]))-3.) <= 1.;
}
VecDoub torusmap(const VecDoub &s) {
VecDoub xx(s);
xx[2] = 0.2*log(5.*s[2]);
return xx;
}

//AV: Function for calculating effective volume of molecules
VecDoub moleculesfunc(const VecDoub &x){
Doub den=1.;
VecDoub f(1);
f[0]=den;
return f;
}
Bool moleculesRegion(const VecDoub &x,int nSites,const MatDoub &Pos){ //pos is array on sig,x,y,z
//x is vector on x,y,z
Bool a=false;
for (int i=0;i<nSites;i++)
if(SQR(x[0]-Pos[i][1])+SQR(x[1]-Pos[i][2])+SQR(x[2]-Pos[i][3])<=SQR(Pos[i][0]))
a=true;
return a;
}


It seems that I have problem in passing arguments to moleculesRegion function. I wonder if you could elaborate on this issue and propose a way to solve this problem. Thanks

I found a solution to this problem by passing the arguments through writing into a txt file and reading from it. Do you have any better solutions for passing arguments through functions?