merge with devel branch

CMakeLists.txt

@@ -11,7 +11,7 @@ FIND_PACKAGE(PkgConfig REQUIRED)
 PKG_CHECK_MODULES(FFTW fftw3 REQUIRED)
 INCLUDE_DIRECTORIES(${FFTW_INCLUDE_DIRS})
 
-FIND_PACKAGE(Boost 1.54 REQUIRED)
+FIND_PACKAGE(Boost 1.43 REQUIRED)
 INCLUDE_DIRECTORIES(${Boost_INCLUDE_DIRS})
 
 FIND_PACKAGE(OpenMP)

bioem.cpp

@@ -57,6 +57,8 @@ int bioem::configure(int ac, char* av[])
 	param.writeAngles = false;
 	param.dumpMap = false;
 	param.loadMap = false;
+	RefMap.readMRC = false;
+	RefMap.readMultMRC = false;
 
 	// *************************************************************************************
 	cout << " ++++++++++++ FROM COMMAND LINE +++++++++++\n\n";
@@ -66,21 +68,26 @@ int bioem::configure(int ac, char* av[])
 
 	try {
 		po::options_description desc("Command line inputs");
-		desc.add_options()
-		("Inputfile", po::value<std::string>(), "Name of input parameter file")
-		("Modelfile", po::value< std::string>() , "Name of model file")
-		("Particlesfile", po::value< std::string>(), "Name of paricles file")
-		("ReadPDB", "(Optional) If reading model file in PDB format")
-		("DumpMaps", "(Optional) Dump maps after they were red from maps file")
-		("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
-		("help", "(Optional) Produce help message")
-		;
+		desc.add_options()
+		("Inputfile", po::value<std::string>(), "(Mandatory) Name of input parameter file")
+		("Modelfile", po::value< std::string>() , "(Mandatory) Name of model file")
+		("Particlesfile", po::value< std::string>(), "(Mandatory) Name of paricles file")
+		("ReadPDB", "(Optional) If reading model file in PDB format")
+		("ReadMRC", "(Optional) If reading particle file in MRC format")
+		("ReadMultipleMRC", "(Optional) If reading Multiple MRCs")
+		("DumpMaps", "(Optional) Dump maps after they were red from maps file")
+		("LoadMapDump", "(Optional) Read Maps from dump instead of maps file")
+		("help", "(Optional) Produce help message")
+		;
+
 
 		po::positional_options_description p;
 		p.add("Inputfile", -1);
 		p.add("Modelfile", -1);
 		p.add("Particlesfile", -1);
 		p.add("ReadPDB", -1);
+		p.add("ReadMRC", -1);
+		p.add("ReadMultipleMRC", -1);
 		p.add("DumpMaps", -1);
 		p.add("LoadMapDump", -1);
 
@@ -118,6 +125,18 @@ int bioem::configure(int ac, char* av[])
 			Model.readPDB = true;
 		}
 
+		if (vm.count("ReadMRC"))
+		{
+			cout << "Reading particle file in MRC format.\n";
+			RefMap.readMRC=true;
+		}
+
+		if (vm.count("ReadMultipleMRC"))
+		{
+			cout << "Reading Multiple MRCs.\n";
+			RefMap.readMultMRC=true;
+		}
+
 		if (vm.count("DumpMaps"))
 		{
 			cout << "Dumping Maps after reading from file.\n";
@@ -360,8 +379,11 @@ int bioem::run()
 
 int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, mycomplex_t* localmultFFT, myfloat_t sumC, myfloat_t sumsquareC, const int startMap)
 {
+	//***************************************************************************************
+	//***** BioEM routine for comparing reference maps to convoluted maps *****
 	if (FFTAlgo)
 	{
+		//With FFT Algorithm
 		#pragma omp parallel for
 		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
 		{
@@ -371,6 +393,7 @@ int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, myc
 	}
 	else
 	{
+		//Without FFT Algorithm
 		#pragma omp parallel for
 		for (int iRefMap = startMap; iRefMap < RefMap.ntotRefMap; iRefMap ++)
 		{
@@ -382,6 +405,9 @@ int bioem::compareRefMaps(int iOrient, int iConv, const myfloat_t* conv_map, myc
 
 inline void bioem::calculateCCFFT(int iRefMap, int iOrient, int iConv, myfloat_t sumC, myfloat_t sumsquareC, mycomplex_t* localConvFFT, mycomplex_t* localCCT, myfloat_t* lCC)
 {
+	//***************************************************************************************
+	//***** Calculating cross correlation in FFTALGOrithm *****
+
 	const mycomplex_t* RefMapFFT = &RefMap.RefMapsFFT[iRefMap * param.FFTMapSize];
 	for(int i = 0; i < param.param_device.NumberPixels * param.param_device.NumberFFTPixels1D; i++)
 	{

include/map.h

@@ -34,8 +34,17 @@ public:
 	int precalculate(bioem_param& param, bioem& bio);
 	int PreCalculateMapsFFT(bioem_param& param);
 
+	int  read_int(int *currlong, FILE *fin, int swap);
+	int  read_float(float *currfloat, FILE *fin, int swap);
+	int  read_float_empty (FILE *fin);
+	int  read_char_float (float *currfloat, FILE *fin) ;
+	int  test_mrc (const char *vol_file, int swap);
+	int  read_MRC(const char* filename,bioem_param& param);
+
 	mycomplex_t* RefMapsFFT;
 
+	bool readMRC,readMultMRC;
+
 	int ntotRefMap;
 	int numPixels;
 	int refMapSize;

map.cpp

@@ -33,9 +33,56 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
 		fread(maps, sizeof(myfloat_t) * refMapSize, ntotRefMap, fp);
 		fclose(fp);
 
-		cout << "Particle Maps read from Map Dump \nTotal Number of particles: " << ntotRefMap ;
-		cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
+		cout << "Particle Maps read from Map Dump\n";
 	}
+	else if(readMRC)
+	{
+		ntotRefMap=0;
+
+		if(readMultMRC)
+		{
+			ifstream input(filemap);
+
+			if (!input.good())
+			{
+				cout << "Failed to open file contaning MRC names: " << filemap << "\n";
+				exit(0);
+			}
+
+			char line[512] = {0};
+			char mapname[100];
+			char tmpm[10] = {0};
+			const char* indifile;
+
+			while (!input.eof())
+			{
+				input.getline(line,512);
+				char tmpVals[100]  = {' '};
+
+				strncpy (tmpVals,line,100);
+				sscanf (tmpVals,"%99c",mapname);
+
+				// Check for last line
+				strncpy (tmpm,mapname,3);		
+				if(strcmp(tmpm,"XXX")!=0)
+				{
+					indifile = mapname;  
+					//Reading Multiple MRC
+					read_MRC(indifile,param);
+
+				}
+				for(int i=0;i<3;i++)mapname[i] = 'X';	
+				for(int i=3;i<100;i++)mapname[i] = {0};	
+
+			}
+			cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";
+			cout << "Particle Maps read from MULTIPLE MRC Files in: " << filemap << "\n" ;
+		} else {
+			read_MRC(filemap,param);
+			cout << "\n++++++++++++++++++++++++++++++++++++++++++ \n";
+			cout << "Particle Maps read from ONE MRC File: " << filemap << "\n" ;
+		}
+	}
 	else
 	{
 		int nummap = -1;
@@ -113,8 +160,7 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
 		cout << ".";
 		ntotRefMap = nummap + 1;
 		maps = (myfloat_t*) reallocchk(maps, refMapSize * sizeof(myfloat_t) * ntotRefMap);
-		cout << "Particle Maps read from Standard File \nTotal Number of particles: " << ntotRefMap ;
-		cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
+		cout << "Particle Maps read from Standard File \n";
 
 		if (param.dumpMap)
 		{
@@ -130,6 +176,9 @@ int bioem_RefMap::readRefMaps(bioem_param& param)
 		}
 	}
 
+	cout << "Total Number of particles: " << ntotRefMap;
+	cout << "\n+++++++++++++++++++++++++++++++++++++++++++ \n";
+
 	sum_RefMap = (myfloat_t*) mallocchk(sizeof(myfloat_t) * ntotRefMap);
 	sumsquare_RefMap = (myfloat_t*) mallocchk(sizeof(myfloat_t) * ntotRefMap);
 
@@ -215,3 +264,267 @@ void bioem_Probability::init(size_t maps, size_t angles, bioem& bio)
 	ptr = bio.malloc_device_host(get_size(maps, angles));
 	set_pointers();
 }
+
+int bioem_RefMap::read_MRC(const char* filename,bioem_param& param)
+{
+	myfloat_t st,st2; 
+	unsigned long count;
+	FILE *fin;
+	float currfloat;
+	int nc, nr, ns, swap, header_ok = 1;
+	float xlen, ylen, zlen;
+	int mode, ncstart, nrstart, nsstart, ispg, nsymbt, lskflg;
+	float a_tmp, b_tmp, g_tmp;
+	int  mx, my, mz,mapc, mapr, maps;
+	float dmin, dmax, dmean;
+	int n_range_viol0, n_range_viol1;
+
+	fin = fopen(filename, "rb");
+	if( fin == NULL ) {
+		cout << "ERROR opening MRC: " << filename;
+		exit(1);
+	}
+	n_range_viol0 = test_mrc(filename,0);
+	n_range_viol1 = test_mrc(filename,1);
+
+
+	if (n_range_viol0 < n_range_viol1) { //* guess endianism 
+		swap = 0;
+		if (n_range_viol0 > 0) {
+			printf(" Warning: %i header field range violations detected in file %s \n", n_range_viol0,filename);
+		}
+	} else {
+		swap = 1;
+		if (n_range_viol1 > 0) {
+			printf("Warning: %i header field range violations detected in file %s \n", n_range_viol1,filename);
+		}
+	}
+	printf("\n+++++++++++++++++++++++++++++++++++++++++++\n");
+	printf("Reading Information from MRC: %s \n", filename);
+	header_ok *= read_int(&nc,fin,swap);
+	header_ok *= read_int(&nr,fin,swap);
+	header_ok *= read_int(&ns,fin,swap);
+	header_ok *= read_int(&mode,fin,swap);
+	header_ok *= read_int(&ncstart,fin,swap);
+	header_ok *= read_int(&nrstart,fin,swap);
+	header_ok *= read_int(&nsstart,fin,swap);
+	header_ok *= read_int(&mx,fin,swap);
+	header_ok *= read_int(&my,fin,swap);
+	header_ok *= read_int(&mz,fin,swap);
+	header_ok *= read_float(&xlen,fin,swap);
+	header_ok *= read_float(&ylen,fin,swap);
+	header_ok *= read_float(&zlen,fin,swap);
+	header_ok *= read_float(&a_tmp,fin,swap);
+	header_ok *= read_float(&b_tmp,fin,swap);
+	header_ok *= read_float(&g_tmp,fin,swap);
+	header_ok *= read_int(&mapc,fin,swap);
+	header_ok *= read_int(&mapr,fin,swap);
+	header_ok *= read_int(&maps,fin,swap);
+	header_ok *= read_float(&dmin,fin,swap);
+	header_ok *= read_float(&dmax,fin,swap);
+	header_ok *= read_float(&dmean,fin,swap);
+	header_ok *= read_int(&ispg,fin,swap);
+	header_ok *= read_int(&nsymbt,fin,swap);
+	header_ok *= read_int(&lskflg,fin,swap);
+
+	printf("Number Columns  = %8d \n",nc);
+	printf("Number Rows     = %8d \n",nr);
+	printf("Number Sections = %8d \n",ns);
+	printf("MODE = %4d (only data type mode 2: 32-bit)\n",mode);
+	printf("NSYMBT = %4d (# bytes symmetry operators)\n",nsymbt);
+
+	/* printf("  NCSTART = %8d  (index of first column, counting from 0)\n",ncstart);
+	printf(">  NRSTART = %8d  (index of first row, counting from 0)\n",nrstart);
+	printf("  NSSTART = %8d  (index of first section, counting from 0)\n",nsstart);
+	printf("       MX = %8d  (# of X intervals in unit cell)\n",mx);
+	printf("       MY = %8d  (# of Y intervals in unit cell)\n",my);
+	printf("       MZ = %8d  (# of Z intervals in unit cell)\n",mz);
+	printf(" X length = %8.3f  (unit cell dimension)\n",xlen);
+	printf(" Y length = %8.3f  (unit cell dimension)\n",ylen);
+	printf(" Z length = %8.3f  (unit cell dimension)\n",zlen);
+	printf("    Alpha = %8.3f  (unit cell angle)\n",a_tmp);
+	printf("     Beta = %8.3f  (unit cell angle)\n",b_tmp);
+	printf("    Gamma = %8.3f  (unit cell angle)\n",g_tmp);
+	printf("     MAPC = %8d  (columns axis: 1=X,2=Y,3=Z)\n",mapc);
+	printf("     MAPR = %8d  (rows axis: 1=X,2=Y,3=Z)\n",mapr);
+	printf("     MAPS = %8d  (sections axis: 1=X,2=Y,3=Z)\n",maps);
+	printf("     DMIN = %8.3f  (minimum density value - ignored)\n",dmin);
+	printf("     DMAX = %8.3f  (maximum density value - ignored)\n",dmax);
+	printf("    DMEAN = %8.3f  (mean density value - ignored)\n",dmean);
+	printf("     ISPG = %8d  (space group number - ignored)\n",ispg);
+	printf("   NSYMBT = %8d  (# bytes storing symmetry operators)\n",nsymbt);
+	printf("   LSKFLG = %8d  (skew matrix flag: 0:none, 1:follows)\n",lskflg);*/
+
+	if (header_ok == 0) {
+		cout << "ERROR reading MRC header: " << filename;
+		exit(1);
+	}
+
+	if(nr!=param.param_device.NumberPixels || nc!=param.param_device.NumberPixels )
+	{
+		cout << "PROBLEM INCONSISTENT NUMBER OF PIXELS IN MAPS AND INPUTFILE ( " << param.param_device.NumberPixels << ", i " << nc << ", j " << nr << ")" << "\n";
+		exit(1);
+	}
+
+	if(mode!=2)
+	{
+		cout << "ERROR with MRC mode " << mode << "\n";
+		cout << "Currently mode 2 is the only one allowed" << "\n";
+		exit(1);
+
+	} else
+	{  
+		rewind (fin);
+		for (count=0; count<256; ++count) if (read_float_empty(fin)==0) {
+			cout << "ERROR Converting Data: " <<  filename;
+			exit(1);
+		}
+
+		for (count=0; count<(unsigned long)nsymbt; ++count) if (read_char_float(&currfloat,fin)==0) {
+			cout << "ERROR Converting Data: " <<  filename;
+			exit(1);    
+		}
+
+		for ( int nmap = 0 ; nmap < ns ; nmap ++ )
+		{
+			st=0.0;
+			st2=0.0;
+			for ( int j = 0 ; j < nr ; j ++ )
+				for ( int i = 0 ; i < nc ; i ++ )
+				{
+					if (read_float(&currfloat,fin,swap)==0) {
+						cout << "ERROR Converting Data: " <<  filename;
+						exit(1);
+					} else {
+						Ref[nmap+ntotRefMap].points[i][j] = currfloat;
+						st += currfloat;
+						st2 += currfloat*currfloat;
+					}
+				}	
+
+				//Normaling maps to zero mean and unit standard deviation
+				st /= float(nr*nc);
+				st2 = sqrt(st2/float(nr*nc)-st*st); 
+				for ( int j = 0 ; j < nr ; j ++ ) for ( int i = 0 ; i < nc ; i ++ ){
+					Ref[nmap+ntotRefMap].points[i][j] =  Ref[nmap+ntotRefMap].points[i][j]/st2 - st/st2;
+					//    if(nmap+ntotRefMap==300)cout << i << " " << j << " " << nmap+ntotRefMap << " " <<  Ref[nmap+ntotRefMap].points[i][j]  << "\n";                          
+				}
+		}
+		ntotRefMap += ns ;
+		//  cout << ntotRefMap << "\n";
+	}
+	fclose (fin);
+
+	return(0);
+}
+
+int bioem_RefMap::read_float(float *currfloat, FILE *fin, int swap) {
+	unsigned char *cptr, tmp;
+
+	if (fread(currfloat,4,1,fin)!=1) return 0;
+	if (swap == 1) {
+		cptr = (unsigned char *)currfloat;
+		tmp = cptr[0];
+		cptr[0]=cptr[3];
+		cptr[3]=tmp;
+		tmp = cptr[1];
+		cptr[1]=cptr[2];
+		cptr[2]=tmp;
+	}
+	return 1;
+}
+
+int  bioem_RefMap::read_int(int *currlong, FILE *fin, int swap) {
+	unsigned char *cptr, tmp;
+
+	if (fread(currlong,4,1,fin)!=1) return 0;
+	if (swap == 1) {
+		cptr = (unsigned char *)currlong;
+		tmp = cptr[0];
+		cptr[0]=cptr[3];
+		cptr[3]=tmp;
+		tmp = cptr[1];
+		cptr[1]=cptr[2];
+		cptr[2]=tmp;
+	}
+	return 1;
+}
+int  bioem_RefMap::read_float_empty (FILE *fin) {
+	float currfloat;
+
+	if (fread(&currfloat,4,1,fin)!=1) return 0;
+	return 1;
+}
+
+int  bioem_RefMap::read_char_float (float *currfloat, FILE *fin) {
+	char currchar;
+
+	if (fread(&currchar,1,1,fin)!=1) return 0;
+	*currfloat=(float)currchar;
+	return 1;
+}
+
+int  bioem_RefMap::test_mrc (const char *vol_file, int swap) {
+	FILE *fin;
+	int nc, nr, ns, mx, my, mz;
+	int mode, ncstart, nrstart, nsstart;
+	float xlen, ylen, zlen;
+	int i, header_ok = 1, n_range_viols = 0;
+	int mapc, mapr, maps;
+	float alpha, beta, gamma;
+	float dmin, dmax, dmean, dummy, xorigin, yorigin, zorigin;
+
+	fin = fopen(vol_file, "rb");
+	if( fin == NULL ) {
+		cout << "ERROR opening MRC: " << vol_file;
+		exit(1);
+	}
+
+	//* read header info 
+	header_ok *= read_int(&nc,fin,swap);
+	header_ok *= read_int(&nr,fin,swap);
+	header_ok *= read_int(&ns,fin,swap);
+	header_ok *= read_int(&mode,fin,swap);
+	header_ok *= read_int(&ncstart,fin,swap);
+	header_ok *= read_int(&nrstart,fin,swap);
+	header_ok *= read_int(&nsstart,fin,swap);
+	header_ok *= read_int(&mx,fin,swap);
+	header_ok *= read_int(&my,fin,swap);
+	header_ok *= read_int(&mz,fin,swap);
+	header_ok *= read_float(&xlen,fin,swap);
+	header_ok *= read_float(&ylen,fin,swap);
+	header_ok *= read_float(&zlen,fin,swap);
+	header_ok *= read_float(&alpha,fin,swap);
+	header_ok *= read_float(&beta,fin,swap);
+	header_ok *= read_float(&gamma,fin,swap);
+	header_ok *= read_int(&mapc,fin,swap);
+	header_ok *= read_int(&mapr,fin,swap);
+	header_ok *= read_int(&maps,fin,swap);
+	header_ok *= read_float(&dmin,fin,swap);
+	header_ok *= read_float(&dmax,fin,swap);
+	header_ok *= read_float(&dmean,fin,swap);
+	for (i=23; i<50; ++i) header_ok *= read_float(&dummy,fin,swap);
+	header_ok *= read_float(&xorigin,fin,swap);
+	header_ok *= read_float(&yorigin,fin,swap);
+	header_ok *= read_float(&zorigin,fin,swap);
+	fclose (fin);
+	if (header_ok == 0) {
+		cout << "ERROR reading MRC header: " << vol_file;
+		exit(1);
+	}
+
+	n_range_viols += (nc>5000); n_range_viols += (nc<0);
+	n_range_viols += (nr>5000); n_range_viols += (nr<0);
+	n_range_viols += (ns>5000); n_range_viols += (ns<0);
+	n_range_viols += (ncstart>5000); n_range_viols += (ncstart<-5000);
+	n_range_viols += (nrstart>5000); n_range_viols += (nrstart<-5000);
+	n_range_viols += (nsstart>5000); n_range_viols += (nsstart<-5000);
+	n_range_viols += (mx>5000); n_range_viols += (mx<0);
+	n_range_viols += (my>5000); n_range_viols += (my<0);
+	n_range_viols += (mz>5000); n_range_viols += (mz<0);
+	n_range_viols += (alpha>360.0f); n_range_viols += (alpha<-360.0f);
+	n_range_viols += (beta>360.0f); n_range_viols += (beta<-360.0f);
+	n_range_viols += (gamma>360.0f); n_range_viols += (gamma<-360.0f);
+
+	return n_range_viols;
+}

model.cpp

@@ -109,8 +109,7 @@ int bioem_model::readModel()
 
 		}
 		nPointsModel = numres;
-		cout << "Protein structure read from PDB\nTotal Number of Residues " << nPointsModel;
-		cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
+		cout << "Protein structure read from PDB\n";
 	}
 	else //Reading model from FILE FORMAT x,y,z,rad,density
 	{
@@ -118,32 +117,35 @@ int bioem_model::readModel()
 		int numres = 0;
 		NormDen = 0.0;
 		FILE *file = fopen ( filemodel , "r" );
-		if ( file != NULL )
+		if (file == NULL)
 		{
-			while ( fgets ( line, sizeof line, file ) != NULL )
+			cout << "Error opening file " << filemode, << "\n";
+			exit(1);
+		}
+		while ( fgets ( line, sizeof line, file ) != NULL )
+		{
+			if (numres >= BIOEM_MODEL_SIZE)
 			{
-				if (numres >= BIOEM_MODEL_SIZE)
-				{
-					cout << "BIOEM_MODEL_SIZE too small\n";
-					exit(1);
-				}
-				float tmpval[5];
-				sscanf(line, "%f %f %f %f %f", &tmpval[0], &tmpval[1], &tmpval[2], &tmpval[3], &tmpval[4]);
-				PointsModel[numres].pos[0] = (myfloat_t) tmpval[0];
-				PointsModel[numres].pos[1] = (myfloat_t) tmpval[1];
-				PointsModel[numres].pos[2] = (myfloat_t) tmpval[2];
-				radiusPointsModel[numres] = (myfloat_t) tmpval[3];
-				densityPointsModel[numres] = (myfloat_t) tmpval[4];
-
-				exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " <<  PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " << densityPointsModel[numres] << "\n";
-				NormDen += densityPointsModel[numres];
-				numres++;
+				cout << "BIOEM_MODEL_SIZE too small\n";
+				exit(1);
 			}
-			nPointsModel = numres;
-			cout << "Protein structure read from Standard File \nTotal Number of Residues " << nPointsModel ;
-			cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
+			float tmpval[5];
+			sscanf(line, "%f %f %f %f %f", &tmpval[0], &tmpval[1], &tmpval[2], &tmpval[3], &tmpval[4]);
+			PointsModel[numres].pos[0] = (myfloat_t) tmpval[0];
+			PointsModel[numres].pos[1] = (myfloat_t) tmpval[1];
+			PointsModel[numres].pos[2] = (myfloat_t) tmpval[2];
+			radiusPointsModel[numres] = (myfloat_t) tmpval[3];
+			densityPointsModel[numres] = (myfloat_t) tmpval[4];
+
+			exampleReadCoor << "RESIDUE " << numres << " " << PointsModel[numres].pos[0] << " " <<  PointsModel[numres].pos[1] << " " << PointsModel[numres].pos[2] << " " << densityPointsModel[numres] << "\n";
+			NormDen += densityPointsModel[numres];
+			numres++;
 		}
+		nPointsModel = numres;
+		cout << "Protein structure read from Standard File\n";
 	}
+	cout << "Total Number of Voxels " << nPointsModel ;
+	cout << "\n+++++++++++++++++++++++++++++++++++++++++ \n";
 	exampleReadCoor.close();
 
 	//Moving to Model to its center of mass:

param.cpp

@@ -42,6 +42,23 @@ int bioem_param::readParameters()
 	// ***************************** Reading Input Parameters ******************************
 	// **************************************************************************************
 
+	// Control for Parameters
+	bool yesPixSi = false;
+	bool yesNumPix = false;
+	bool yesGPal = false;
+	bool yesGPbe = false;
+	bool yesGPEnv = false;
+	bool yesGPamp = false;
+	bool yesGPpha = false;
+	bool yesSTEnv = false;
+	bool yesSTamp = false;
+	bool yesSTpha = false;
+	bool yesGSPamp = false ;
+	bool yesGSPEnv = false ;
+	bool yesGSPpha = false ;
+	bool yesMDC = false ;
+	bool yesGCen = false ;	
+
 	ifstream input(fileinput);
 	if (!input.good())
 	{
@@ -70,116 +87,162 @@ int bioem_param::readParameters()
 		{
 			token = strtok(NULL, " ");
 			pixelSize = atof(token);
-			cout << "Pixel Sixe " << pixelSize << "\n";
-
+			if (pixelSize < 0 ) { cout << "*** Error: Negative pixelSize "; exit(1);}
+			cout << "Pixel Sixe " << pixelSize << "\n";
+			yesPixSi= true;
 		}
 		else if (strcmp(token, "NUMBER_PIXELS") == 0)
 		{
 			token = strtok(NULL, " ");
 			param_device.NumberPixels = int(atoi(token));
-			cout << "Number of Pixels " << param_device.NumberPixels << "\n";
-
+			if (param_device.NumberPixels < 0 ) { cout << "*** Error: Negative Number of Pixels "; exit(1);}
+			cout << "Number of Pixels " << param_device.NumberPixels << "\n";			
+			yesNumPix= true ;
 		}
 		else if (strcmp(token, "GRIDPOINTS_ALPHA") == 0)
 		{
 			token = strtok(NULL, " ");
 			angleGridPointsAlpha = int(atoi(token));
-			cout << "Grid points alpha " << angleGridPointsAlpha << "\n";
-
+			if (angleGridPointsAlpha < 0 ) { cout << "*** Error: Negative GRIDPOINTS_ALPHA "; exit(1);}
+			cout << "Grid points alpha " << angleGridPointsAlpha << "\n";
+			yesGPal= true;
 		}
 		else if (strcmp(token, "GRIDPOINTS_BETA") == 0)
 		{
 			token = strtok(NULL, " ");
 			angleGridPointsBeta = int(atoi(token));
-			cout << "Grid points beta " << angleGridPointsBeta << "\n";
-
+			if (angleGridPointsBeta < 0 ) { cout << "*** Error: Negative GRIDPOINTS_BETA "; exit(1);}
+			cout << "Grid points beta " << angleGridPointsBeta << "\n";
+			yesGPbe= true;
 		}
 
 		else if (strcmp(token, "GRIDPOINTS_ENVELOPE") == 0)
 		{
 			token = strtok(NULL, " ");
 			numberGridPointsEnvelop = int(atoi(token));
-			cout << "Grid points envelope " << numberGridPointsEnvelop << "\n";
-
+			if (numberGridPointsDisplaceCenter < 0 ) { cout << "*** Error: Negative GRIDPOINTS_ENVELOPE "; exit(1);}
+			cout << "Grid points envelope " << numberGridPointsEnvelop << "\n";
+			yesGPEnv = true;
 		}
 		else if (strcmp(token, "START_ENVELOPE") == 0)
 		{
 			token = strtok(NULL, " ");
 			startGridEnvelop = atof(token);
-			cout << "Start Envelope " << startGridEnvelop << "\n";
-