Configuration.cpp

#include "Configuration.h"
#include "Angle.h"
#include "Dihedral.h"
#include "Restraint.h"
#include <cmath>
#include <cassert>
#include <stdlib.h>     /* srand, rand */
#include <time.h>       /* time */
#include <gsl/gsl_rng.h>
#include <gsl/gsl_randist.h>
#include <gsl/gsl_math.h>

#define gpuErrchk(ans) { gpuAssert((ans), __FILE__, __LINE__); }
inline void gpuAssert(cudaError_t code, const char *file, int line, bool abort=true) {
   if (code != cudaSuccess) {
      fprintf(stderr,"CUDA Error: %s %s %d\n", cudaGetErrorString(code), __FILE__, line);
      if (abort) exit(code);
   }
}

Configuration::Configuration(const char* config_file, int simNum, bool debug) :
		simNum(simNum) {
	// Read the parameters.
	//type_d = NULL;
	kTGrid_d = NULL;
	//bonds_d = NULL;
	//bondMap_d = NULL;
	//excludes_d = NULL;
	//excludeMap_d = NULL;
	//angles_d = NULL;
	//dihedrals_d = NULL;
	setDefaults();
	readParameters(config_file);

	// Get the number of particles
	// printf("\nCounting particles specified in the ");
	if (restartCoordinates.length() > 0) {
    // Read them from the restart file.
	    // printf("restart file.\n");
		num = countRestart(restartCoordinates.val());
  } else {
    if (readPartsFromFile) readAtoms();
    if (numPartsFromFile > 0) {
      // Determine number of particles from input file (PDB-style)
	// printf("input file.\n");
      num = numPartsFromFile;
    } else {
      // Sum up all particles in config file
	// printf("configuration file.\n");
      //int num0 = 0;
      num = 0;
      for (int i = 0; i < numParts; i++) num += part[i].num;
      //num = num0;
    }
  } // end result: variable "num" is set


	// Set the number capacity
	printf("\n%d particles\n", num);
	if (numCap <= 0) numCap = numCapFactor*num; // max number of particles
	if (numCap <= 0) numCap = 20;

	// Allocate particle variables.
	pos = new Vector3[num * simNum];
        //Han-Yi Chou
        if (ParticleDynamicType == String("Langevin") || ParticleDynamicType == String("NoseHooverLangevin"))
            momentum = new Vector3[num * simNum];

	type   = new int[num * simNum];
	serial = new int[num * simNum];
	posLast = new Vector3[num * simNum];
        //Han-Yi Chou
        if(ParticleDynamicType == String("Langevin") || ParticleDynamicType == String("NoseHooverLangevin"))
           momLast = new Vector3[num * simNum];
	name = new String[num * simNum];
	currSerial = 0;


  // Now, load the coordinates
	loadedCoordinates = false;
        loadedMomentum    = false; //Han-Yi Chou

  //I need kT here Han-Yi Chou
  kT = temperature * 0.0019872065f; // `units "k K" "kcal_mol"`
  //kT = temperature * 0.593f;
 // If we have a restart file - use it
	if (restartCoordinates.length() > 0) {
		loadRestart(restartCoordinates.val()); 
		printf("Loaded %d restart coordinates from `%s'.\n", num, restartCoordinates.val());
		printf("Particle numbers specified in the configuration file will be ignored.\n");
		loadedCoordinates = true;
                //Han-Yi Chou Langevin dynamic
                if (ParticleDynamicType == String("Langevin") || ParticleDynamicType == String("NoseHooverLangevin"))
                {
                    if (restartMomentum.length() > 0)
                    {
                        loadRestartMomentum(restartMomentum.val());
                        printf("Loaded %d restart momentum from `%s'.\n", num, restartMomentum.val());
                        printf("Particle numbers specified in the configuration file will be ignored.\n");
                        loadedMomentum = true;
                    }
                    else
                    {
                        printf("Warning: There is no restart momentum file when using restart coordinates in Langevin Dynamics\n");
                        printf("Initialize with Boltzmann distribution\n");
                        loadedMomentum = Boltzmann(COM_Velocity, num * simNum);
                    }
               }
	} 
        else 
        {
		// Load coordinates from a file?
		if (numPartsFromFile > 0) {
			loadedCoordinates = true;
			for (int i = 0; i < num; i++) {
				int numTokens = partsFromFile[i].tokenCount();

				// Break the line down into pieces (tokens) so we can process them individually
				String* tokenList = new String[numTokens];
				partsFromFile[i].tokenize(tokenList);

				int currType = 0;
				for (int j = 0; j < numParts; j++)
					if (tokenList[2] == part[j].name)
						currType = j;

				for (int s = 0; s < simNum; ++s)
				    type[i + s*num] = currType;

				serial[i] = currSerial++;

				pos[i] = Vector3(atof(tokenList[3].val()), atof(tokenList[4].val()), atof(tokenList[5].val()));
                                //Han-Yi Chou
                                if (ParticleDynamicType == String("Langevin") || ParticleDynamicType == String("NoseHooverLangevin"))
                                {
                                    loadedMomentum = true;
                                    if(numTokens == 9)
                                        momentum[i] = Vector3(atof(tokenList[6].val()), atof(tokenList[7].val()), atof(tokenList[8].val()));
                                    else
                                    {
                                        printf("Error occurs in %s at line %d. Please specify momentum\n", __FILE__, __LINE__);
                                        assert(1==2);
                                    }
                                }
			}
			delete[] partsFromFile;
			partsFromFile = NULL;
                        //Han-Yi Chou
                        for(int i = 1; i < simNum; ++i)
                            for(int j = 0; j < num; ++j)
                                serial[j + num * i] = currSerial++;
                }
                else 
                {
	            // Not loading coordinates from a file
	            populate();
	            if (inputCoordinates.length() > 0) 
                    {
		        printf("Loading coordinates from %s ... ", inputCoordinates.val());
			loadedCoordinates = loadCoordinates(inputCoordinates.val());
			if (loadedCoordinates)
			    printf("done!\n");
	            }
                    if(ParticleDynamicType == String("Langevin") || ParticleDynamicType == String("NoseHooverLangevin"))
                    {
                        if (inputMomentum.length() > 0) 
                        {
                            printf("Loading momentum from %s ... ", inputMomentum.val());
                            loadedMomentum = loadMomentum(inputMomentum.val());
                            if (loadedMomentum)
                                printf("done!\n");
                        }
                        else
                            loadedMomentum = Boltzmann(COM_Velocity, num * simNum);
                    }
                }
            }
        //Check initialize momentum
        //if(ParticleDynamicType == String("Langevin"))
            //PrintMomentum();
	/* Initialize exclusions */
	excludeCapacity = 256;
	numExcludes = 0;
	excludes = new Exclude[excludeCapacity];

	if (readExcludesFromFile) readExcludes();
	if (readBondsFromFile) readBonds();
	if (readAnglesFromFile) readAngles();
	if (readDihedralsFromFile) readDihedrals();
	if (readRestraintsFromFile) readRestraints();

	if (temperatureGridFile.length() != 0) {
		printf("\nFound temperature grid file: %s\n", temperatureGridFile.val());
		tGrid = new BaseGrid(temperatureGridFile.val());
		printf("Loaded `%s'.\n", temperatureGridFile.val());
		printf("Grid size %s.\n", tGrid->getExtent().toString().val());

		// TODO: ask Max Belkin what this is about and how to remove hard-coded temps
		float ToSo = 1.0f / (295.0f * 4.634248239f); // 1 / (To * sigma(To))
		sigmaT = new BaseGrid(*tGrid);
		sigmaT->shift(-122.8305f);
		sigmaT->scale(0.0269167f);
		sigmaT->mult(*tGrid);
		sigmaT->scale(ToSo);

		kTGrid = new BaseGrid(*tGrid);
		float factor = 0.0019872065f; // `units "k K" "kcal_mol"`
		kTGrid->scale(factor);
		// char outFile[256];
		// char comment[256]; sprintf(comment,"KTGrid");
		// sprintf(outFile,"kTGrid.dx");
		// kTGrid->write(outFile, comment);
	}

	printf("\nFound %d particle types.\n", numParts);
	// Load the potential grids.
	printf("Loading the potential grids...\n");
	for (int i = 0; i < numParts; i++) {
		// Decide which type of grid is given.
		String map = partGridFile[i];
		int len = map.length();
		if (len >= 3 && map[len-3]=='.' && map[len-2]=='d' && map[len-1]=='x') {
			// A dx file. Load the old-fashioned way.
			part[i].pmf = new BaseGrid(map.val());
			if (partGridFileScale[i] != 1.0f) part[i].pmf->scale(partGridFileScale[i]);

			part[i].meanPmf = part[i].pmf->mean();
			printf("Loaded dx potential grid `%s'.\n", map.val());
			printf("Grid size %s.\n", part[i].pmf->getExtent().toString().val());
		} else if  (len >= 4 && map[len-4]=='.' && map[len-3]=='d' && map[len-2]=='e' && map[len-1]=='f') {
			// A system definition file.
			String rootGrid = OverlordGrid::readDefFirst(map);
			OverlordGrid* over = new OverlordGrid(rootGrid.val());
			int count = over->readDef(map);
			printf("Loaded system def file `%s'.\n", map.val());
			printf("Found %d unique grids.\n", over->getUniqueGridNum());
			printf("Linked %d subgrids.\n", count);

			part[i].pmf = static_cast<BaseGrid*>(over);
			part[i].meanPmf = part[i].pmf->mean();
		} else {
			printf("WARNING: Unrecognized gridFile extension. Must be *.def or *.dx.\n");
			exit(-1);
		}

		if (partForceXGridFile[i].length() != 0) {
			part[i].forceXGrid = new BaseGrid(partForceXGridFile[i].val());
			printf("Loaded `%s'.\n", partForceXGridFile[i].val());
			printf("Grid size %s.\n", part[i].forceXGrid->getExtent().toString().val());
		}

		if (partForceYGridFile[i].length() != 0) {
			part[i].forceYGrid = new BaseGrid(partForceYGridFile[i].val());
			printf("Loaded `%s'.\n", partForceYGridFile[i].val());
			printf("Grid size %s.\n", part[i].forceYGrid->getExtent().toString().val());
		}

		if (partForceZGridFile[i].length() != 0) {
			part[i].forceZGrid = new BaseGrid(partForceZGridFile[i].val());
			printf("Loaded `%s'.\n", partForceZGridFile[i].val());
			printf("Grid size %s.\n", part[i].forceZGrid->getExtent().toString().val());
		}

		if (partDiffusionGridFile[i].length() != 0) {
			part[i].diffusionGrid = new BaseGrid(partDiffusionGridFile[i].val());
			printf("Loaded `%s'.\n", partDiffusionGridFile[i].val());
			printf("Grid size %s.\n", part[i].diffusionGrid->getExtent().toString().val());
		}

		if (temperatureGridFile.length() != 0) {
			if (partDiffusionGridFile[i].length() != 0) {
				part[i].diffusionGrid->mult(*sigmaT);
			} else {
				part[i].diffusionGrid = new BaseGrid(*sigmaT);
				part[i].diffusionGrid->scale(part[i].diffusion);
				// char outFile[256];
				// char comment[256]; sprintf(comment,"Diffusion for particle type %d", i);
				// sprintf(outFile,"diffusion%d.dx",i);
				// part[i].diffusionGrid->write(outFile, comment);
			}
		}
	}

    // Load reservoir files if any
    for (int i = 0; i < numParts; i++) {
        if (partReservoirFile[i].length() != 0) {
            printf("\nLoading the reservoirs for %s... \n", part[i].name.val());
            part[i].reservoir = new Reservoir(partReservoirFile[i].val());
            int nRes = part[i].reservoir->length();
            printf("\t -> %d reservoir(s) found in `%s'.\n", nRes, partReservoirFile[i].val());
        }
    }

    // Get the system dimensions
    // from the dimensions of supplied 3D potential maps
    if (size.length2() > 0) {	// use size if it's defined
	if (basis1.length2() > 0 || basis2.length2() > 0 || basis3.length2() > 0)
	    printf("WARNING: both 'size' and 'basis' were specified... using 'size'\n"); 
	basis1 = Vector3(size.x,0,0);
	basis2 = Vector3(0,size.y,0);
	basis3 = Vector3(0,0,size.z);
    }
    if (basis1.length2() > 0 && basis2.length2() > 0 && basis3.length2() > 0) {
	sys = new BaseGrid( Matrix3(basis1,basis2,basis3), origin, 1, 1, 1 );
    } else {
	// TODO: use largest system in x,y,z
	sys = part[0].pmf;
    }
    sysDim = sys->getExtent();

// RBTODO: clean this mess up
	/* // RigidBodies... */
	/* if (numRigidTypes > 0) { */
	/* 	printf("\nCounting rigid bodies specified in the configuration file.\n"); */
	/* 	numRB = 0; */

	/* 	// grow list of rbs */
	/* 	for (int i = 0; i < numRigidTypes; i++) {			 */
	/* 		numRB += rigidBody[i].num; */

	/* 		std::vector<RigidBody> tmp; */
	/* 		for (int j = 0; j < rigidBody[i].num; j++) { */
	/* 			tmp.push_back( new RigidBody( this, rigidBody[i] ) ); */
	/* 		} */

	/* 		rbs.push_back(tmp); */
	/* 	} */
		// // state data
		// rbPos = new Vector3[numRB * simNum];
		// type = new int[numRB * simNum];

	/* } */
	/* printf("Initial RigidBodies: %d\n", numRB); */


	// Create exclusions from the exclude rule, if it was specified in the config file
	if (excludeRule != String("")) {
		int oldNumExcludes = numExcludes;
		Exclude* newExcludes = makeExcludes(bonds, bondMap, num, numBonds, excludeRule, numExcludes);
		if (excludes == NULL) {
			excludes = new Exclude[numExcludes];
		} else if (numExcludes >= excludeCapacity) {
			Exclude* tempExcludes = excludes;
			excludes = new Exclude[numExcludes];
			for (int i = 0; i < oldNumExcludes; i++)
				excludes[i] = tempExcludes[i];
			delete tempExcludes;
		}
		for (int i = oldNumExcludes; i < numExcludes; i++)
			excludes[i] = newExcludes[i - oldNumExcludes];
		printf("Built %d exclusions.\n",numExcludes);
	}

	buildExcludeMap();

	// Count number of particles of each type
	numPartsOfType = new int[numParts];
	for (int i = 0; i < numParts; ++i) {
		numPartsOfType[i] = 0;
	}
	for (int i = 0; i < num; ++i) {
		++numPartsOfType[type[i]];
	}

	// Some geometric stuff that should be gotten rid of.
	Vector3 buffer = (sys->getCenter() + 2.0f*sys->getOrigin())/3.0f;
	initialZ = buffer.z;

	// Set the initial conditions.
	// Do the initial conditions come from restart coordinates?
	// inputCoordinates are ignored if restartCoordinates exist.
	/*
	if (restartCoordinates.length() > 0) {
		loadRestart(restartCoordinates.val());
		printf("Loaded %d restart coordinates from `%s'.\n", num, restartCoordinates.val());
		printf("Particle numbers specified in the configuration file will be ignored.\n");
	} else {
		// Set the particle types.

		// Load coordinates from a file?
		if (numPartsFromFile > 0) {
			for (int i = 0; i < num; i++) {
				int numTokens = partsFromFile[i].tokenCount();

				// Break the line down into pieces (tokens) so we can process them individually
				String* tokenList = new String[numTokens];
				partsFromFile[i].tokenize(tokenList);
				int currType = 0;
				for (int j = 0; j < numParts; j++)
					if (tokenList[2] == part[j].name)
						currType = j;
				type[i] = currType;
				serial[i] = currSerial;
				currSerial++;

				pos[i] = Vector3(atof(tokenList[3].val()), atof(tokenList[4].val()), atof(tokenList[5].val()));
			}
			if (partsFromFile != NULL) {
				delete[] partsFromFile;
				partsFromFile = NULL;
			}
		} else if (inputCoordinates.length() > 0) {
			populate();
			printf("Loading coordinates from %s.\n", inputCoordinates.val());
			bool loaded = loadCoordinates(inputCoordinates.val());
			if (loaded) 
				printf("Loaded initial coordinates from %s.\n", inputCoordinates.val());
		}
	}
	*/
	

	// Get the maximum particle radius.
	minimumSep = 0.0f;
	for (int i = 0; i < numParts; ++i)
		minimumSep = std::max(minimumSep, part[i].radius);
	minimumSep *= 2.5f; // Make it a little bigger.

	// Default outputEnergyPeriod
	if (outputEnergyPeriod < 0)
		outputEnergyPeriod = 10 * outputPeriod;
	
	// If we are running with debug ON, ask the user which force computation to use
	if (debug)
		getDebugForce();

	printf("\n");
	switchStart = cutoff - switchLen;
	if (fullLongRange == 0)
		printf("Cutting off the potential from %.10g to %.10g.\n", switchStart, switchStart+switchLen);
	
	if (fullLongRange != 0)
		printf("No cell decomposition created.\n");

}

Configuration::~Configuration() {
	// System state
	delete[] pos;
        //Han-Yi Chou
        if (ParticleDynamicType == String("Langevin") || ParticleDynamicType == String("NoseHooverLangevin"))
            delete[] momentum;

	delete[] posLast;
        //Han-Yi Chou
        if (ParticleDynamicType == String("Langevin") || ParticleDynamicType == String("NoseHooverLangevin"))
            delete[] momLast;

	delete[] type;
	delete[] name;
	
	// Particle parameters
	delete[] part;
	delete[] partGridFile;
	delete[] partGridFileScale;
	delete[] partForceXGridFile;
	delete[] partForceYGridFile;
	delete[] partForceZGridFile;
	delete[] partDiffusionGridFile;
	delete[] partReservoirFile;
	partRigidBodyGrid.clear();
	
	// TODO: plug memory leaks
	if (partsFromFile != NULL) delete[] partsFromFile;
	if (bonds != NULL) delete[] bonds;
	if (bondMap != NULL) delete[] bondMap;
	if (excludes != NULL) delete[] excludes;
	if (excludeMap != NULL) delete[] excludeMap;
	if (angles != NULL) delete[] angles;
	if (dihedrals != NULL) delete[] dihedrals;

	delete[] numPartsOfType;
	  
	// Table parameters
	delete[] partTableFile;
	delete[] partTableIndex0;
	delete[] partTableIndex1;

	delete[] bondTableFile;

	delete[] angleTableFile;

	delete[] dihedralTableFile;

	//if (type_d != NULL) {
		//gpuErrchk(cudaFree(type_d));
		gpuErrchk(cudaFree(sys_d));
		gpuErrchk(cudaFree(kTGrid_d));
		gpuErrchk(cudaFree(part_d));
		//gpuErrchk(cudaFree(bonds_d));
		//gpuErrchk(cudaFree(bondMap_d));
		//gpuErrchk(cudaFree(excludes_d));
		//gpuErrchk(cudaFree(excludeMap_d));
		//gpuErrchk(cudaFree(angles_d));
		//gpuErrchk(cudaFree(dihedrals_d));
	//}
}

void Configuration::copyToCUDA() {
	printf("Copying particle data to GPU %d\n", GPUManager::current());

	BrownianParticleType **part_addr = new BrownianParticleType*[numParts];

	// Copy the BaseGrid objects and their member variables/objects
	gpuErrchk(cudaMalloc(&part_d, sizeof(BrownianParticleType*) * numParts));
	// TODO: The above line fails when there is not enough memory. If it fails, stop.
	
	for (int i = 0; i < numParts; i++) {
		BaseGrid *pmf = NULL, *diffusionGrid = NULL;
		BrownianParticleType *b = new BrownianParticleType(part[i]);
		// Copy PMF
		if (part[i].pmf != NULL) {
			float *val = NULL;
			size_t sz = sizeof(float) * part[i].pmf->getSize();
		  gpuErrchk(cudaMalloc(&pmf, sizeof(BaseGrid)));
		  gpuErrchk(cudaMalloc(&val, sz));
		  gpuErrchk(cudaMemcpyAsync(val, part[i].pmf->val, sz, cudaMemcpyHostToDevice));
		  BaseGrid *pmf_h = new BaseGrid(*part[i].pmf);
			pmf_h->val = val;
			gpuErrchk(cudaMemcpy(pmf, pmf_h, sizeof(BaseGrid), cudaMemcpyHostToDevice));
			pmf_h->val = NULL;
		}
		
		// Copy the diffusion grid
		if (part[i].diffusionGrid != NULL) {
			float *val = NULL;
			size_t sz = sizeof(float) * part[i].diffusionGrid->getSize();
		  BaseGrid *diffusionGrid_h = new BaseGrid(*part[i].diffusionGrid);
		  gpuErrchk(cudaMalloc(&diffusionGrid, sizeof(BaseGrid)));
		  gpuErrchk(cudaMalloc(&val, sz));
			diffusionGrid_h->val = val;
			gpuErrchk(cudaMemcpyAsync(diffusionGrid, diffusionGrid_h, sizeof(BaseGrid),
					cudaMemcpyHostToDevice));
		  gpuErrchk(cudaMemcpy(val, part[i].diffusionGrid->val, sz, cudaMemcpyHostToDevice));
			diffusionGrid_h->val = NULL;
		}
		
		b->pmf = pmf;
		b->diffusionGrid = diffusionGrid;
		gpuErrchk(cudaMalloc(&part_addr[i], sizeof(BrownianParticleType)));
		gpuErrchk(cudaMemcpyAsync(part_addr[i], b, sizeof(BrownianParticleType),
				cudaMemcpyHostToDevice));
	}
	// RBTODO: moved this out of preceding loop; was that correct?
	gpuErrchk(cudaMemcpyAsync(part_d, part_addr, sizeof(BrownianParticleType*) * numParts,
				cudaMemcpyHostToDevice));

	// kTGrid_d
	kTGrid_d = NULL;
	if (temperatureGridFile.length() > 0) {
		gpuErrchk(cudaMalloc(&kTGrid_d, sizeof(BaseGrid)));
		gpuErrchk(cudaMemcpyAsync(kTGrid_d, kTGrid, sizeof(BaseGrid), cudaMemcpyHostToDevice));
	}

	// type_d and sys_d
	gpuErrchk(cudaMalloc(&sys_d, sizeof(BaseGrid)));
	gpuErrchk(cudaMemcpyAsync(sys_d, sys, sizeof(BaseGrid), cudaMemcpyHostToDevice));
	/*gpuErrchk(cudaMalloc(&type_d, sizeof(int) * num * simNum));
	gpuErrchk(cudaMemcpyAsync(type_d, type, sizeof(int) * num * simNum, cudaMemcpyHostToDevice));
	
	if (numBonds > 0) {
		// bonds_d
		gpuErrchk(cudaMalloc(&bonds_d, sizeof(Bond) * numBonds));
		gpuErrchk(cudaMemcpyAsync(bonds_d, bonds, sizeof(Bond) * numBonds, cudaMemcpyHostToDevice));
		
		// bondMap_d
		gpuErrchk(cudaMalloc(&bondMap_d, sizeof(int2) * num));
		gpuErrchk(cudaMemcpyAsync(bondMap_d, bondMap, sizeof(int2) * num, cudaMemcpyHostToDevice));
	}

	if (numExcludes > 0) {
		// excludes_d
		gpuErrchk(cudaMalloc(&excludes_d, sizeof(Exclude) * numExcludes));
		gpuErrchk(cudaMemcpyAsync(excludes_d, excludes, sizeof(Exclude) * numExcludes,
				cudaMemcpyHostToDevice));
		
		// excludeMap_d
		gpuErrchk(cudaMalloc(&excludeMap_d, sizeof(int2) * num));
		gpuErrchk(cudaMemcpyAsync(excludeMap_d, excludeMap, sizeof(int2) * num,
				cudaMemcpyHostToDevice));
	}

	if (numAngles > 0) {
		// angles_d
		gpuErrchk(cudaMalloc(&angles_d, sizeof(Angle) * numAngles));
		gpuErrchk(cudaMemcpyAsync(angles_d, angles, sizeof(Angle) * numAngles,
				cudaMemcpyHostToDevice));
	}

	if (numDihedrals > 0) {
		// dihedrals_d
		gpuErrchk(cudaMalloc(&dihedrals_d, sizeof(Dihedral) * numDihedrals));
		gpuErrchk(cudaMemcpyAsync(dihedrals_d, dihedrals,
												 		  sizeof(Dihedral) * numDihedrals,
														 	cudaMemcpyHostToDevice));
	}*/
	gpuErrchk(cudaDeviceSynchronize());
}

void Configuration::setDefaults() {
    // System parameters
	outputName = "out";
	timestep = 1e-5f;
	rigidBodyGridGridPeriod = 1;
	steps = 100;
	seed = 0;

	origin = Vector3(0,0,0);
	size = Vector3(0,0,0);
	basis1 = Vector3(0,0,0);
	basis2 = Vector3(0,0,0);
	basis3 = Vector3(0,0,0);
	
	inputCoordinates = "";
	restartCoordinates = "";
        //Han-Yi Chou
        inputMomentum = "";
        restartMomentum = "";
	copyReplicaCoordinates = 1;
	numberFluct = 0;
	numberFluctPeriod = 200;
	interparticleForce = 1;
	tabulatedPotential = 0;
	fullLongRange = 1;
	//	kTGridFile = ""; // Commented out for an unknown reason
	temperature = 295.0f;
	temperatureGridFile = "";
	coulombConst = 566.440698f/92.0f;
	electricField = 0.0f;
	cutoff = 10.0f;
	switchLen = 2.0f;
	pairlistDistance = 2.0f;
	imdForceScale = 1.0f;
	outputPeriod = 200;
	outputEnergyPeriod = -1;
	outputFormat = TrajectoryWriter::formatDcd;
	currentSegmentZ = -1.0f;
	numCap = 0;
	decompPeriod = 10;
	readPartsFromFile = 0;
	numPartsFromFile = 0;
	partsFromFile = NULL;
	readBondsFromFile = false;
	numBonds = 0;
	bonds = NULL;
	bondMap = NULL;
	numTabBondFiles = 0;
	readExcludesFromFile = false;
	numExcludes = 0;
	excludeCapacity = 256;
	excludes = NULL;
	excludeMap = NULL;
	excludeRule = "";
	readAnglesFromFile = false;
	numAngles = 0;
	angles = NULL;
	numTabAngleFiles = 0;
	readDihedralsFromFile = false;
	numDihedrals = 0;
	dihedrals = NULL;
	numTabDihedralFiles = 0;

	readRestraintsFromFile = false;
	numRestraints = 0;
	restraints = NULL;

        //Han-Yi Chou default values
        ParticleDynamicType  = String("Brown");
        RigidBodyDynamicType = String("Brown");
        COM_Velocity = Vector3(0.f,0.f,0.f);
        ParticleLangevinIntegrator = String("BAOAB"); //The default is BAOAB

	// Hidden parameters
	// Might be parameters later
	numCapFactor = 5;
}

int Configuration::readParameters(const char * config_file) {
	Reader config(config_file);
	printf("Read config file %s\n", config_file);

	// Get the number of particles.
	const int numParams = config.length();
	numParts = config.countParameter("particle");
	numRigidTypes = config.countParameter("rigidBody");

	// Allocate the particle variables.
	part = new BrownianParticleType[numParts];
	partGridFile = new String[numParts];
	partGridFileScale = new float[numParts];
	partForceXGridFile = new String[numParts];
	partForceYGridFile = new String[numParts];
	partForceZGridFile = new String[numParts];
	partDiffusionGridFile = new String[numParts];
	partReservoirFile = new String[numParts];
	partRigidBodyGrid.resize(numParts);
	
	// Allocate the table variables.
	partTableFile = new String[numParts*numParts];
	partTableIndex0 = new int[numParts*numParts];
	partTableIndex1 = new int[numParts*numParts];

	// Allocate rigid body types
	rigidBody = new RigidBodyType[numRigidTypes];
	
	// Set a default
	for (int i = 0; i < numParts; ++i) {
	    partGridFileScale[i] = 1.0f;
	}

	int btfcap = 10;
	bondTableFile = new String[btfcap];

	int atfcap = 10;
	angleTableFile = new String[atfcap];

	int dtfcap = 10;
	dihedralTableFile = new String[dtfcap];

	int currPart = -1;
	int currTab = -1;
	int currBond = -1;
	int currAngle = -1;
	int currDihedral = -1;
	int currRB = -1;

	int partClassPart =  0;
	int partClassRB   =  1;
	int currPartClass = -1;				// 0 => particle, 1 => rigidBody
	
	for (int i = 0; i < numParams; i++) {
		String param = config.getParameter(i);
		String value = config.getValue(i);
		// GLOBAL
		if (param == String("outputName"))
			outputName = value;
		else if (param == String("timestep"))
			timestep = (float) strtod(value.val(), NULL);
		else if (param == String("rigidBodyGridGridPeriod"))
			rigidBodyGridGridPeriod = atoi(value.val());
		else if (param == String("steps"))
			steps = atol(value.val());
		else if (param == String("seed"))
			seed = atoi(value.val());
		else if (param == String("origin"))
		    origin = stringToVector3( value );
		else if (param == String("systemSize"))
		    size = stringToVector3( value );
		else if (param == String("basis1"))
		    basis1 = stringToVector3( value );
		else if (param == String("basis2"))
		    basis2 = stringToVector3( value );
		else if (param == String("basis3"))
		    basis3 = stringToVector3( value );
		else if (param == String("inputCoordinates"))
			inputCoordinates = value;
		else if (param == String("restartCoordinates"))
			restartCoordinates = value;

                //Han-Yi Chou
                else if (param == String("inputMomentum"))
                        inputMomentum = value;
                else if (param == String("restartMomentum"))
                        restartMomentum = value;
		else if (param == String("copyReplicaCoordinates"))
		        copyReplicaCoordinates = atoi(value.val());
		else if (param == String("temperature"))
			temperature =  (float) strtod(value.val(),NULL);
		else if (param == String("temperatureGrid"))
			temperatureGridFile = value;
		else if (param == String("numberFluct"))
			numberFluct = atoi(value.val());
		else if (param == String("numberFluctPeriod"))
			numberFluctPeriod = atoi(value.val());
		else if (param == String("interparticleForce"))
			interparticleForce = atoi(value.val());
		else if (param == String("fullLongRange") || param == String("fullElect") )
			fullLongRange = atoi(value.val());
		else if (param == String("coulombConst"))
			coulombConst = (float) strtod(value.val(), NULL);
		else if (param == String("electricField"))
			electricField = (float) strtod(value.val(), NULL);
		else if (param == String("cutoff"))
			cutoff = (float) strtod(value.val(), NULL);
		else if (param == String("switchLen"))
			switchLen = (float) strtod(value.val(), NULL);
		else if (param == String("pairlistDistance"))
			pairlistDistance = (float) strtod(value.val(), NULL);
		else if (param == String("scaleIMDForce"))
			imdForceScale = (float) strtod(value.val(), NULL);		
		else if (param == String("outputPeriod"))
			outputPeriod = atoi(value.val());
		else if (param == String("outputEnergyPeriod"))
			outputEnergyPeriod = atoi(value.val());
		else if (param == String("outputFormat"))
			outputFormat = TrajectoryWriter::getFormatCode(value);
		else if (param == String("currentSegmentZ"))
			currentSegmentZ = (float) strtod(value.val(), NULL);
		else if (param == String("numCap"))
			numCap = atoi(value.val());
		else if (param == String("decompPeriod"))
			decompPeriod = atoi(value.val());

                //Han-Yi Chou
                else if (param == String("ParticleDynamicType"))
                    ParticleDynamicType = value;
                else if (param == String("RigidBodyDynamicType"))
                    RigidBodyDynamicType = value;
                else if (param == String("ParticleLangevinIntegrator"))
                    ParticleLangevinIntegrator = value;

		// PARTICLES
		else if (param == String("particle")) {
			part[++currPart] = BrownianParticleType(value);
			currPartClass = partClassPart;
		}
                else if (param == String("mu")) // for Nose-Hoover Langevin
                        part[currPart].mu = (float) strtod(value.val(), NULL);
		else if (param == String("forceXGridFile"))
			partForceXGridFile[currPart] = value;
		else if (param == String("forceYGridFile"))
			partForceYGridFile[currPart] = value;
		else if (param == String("forceZGridFile"))
			partForceZGridFile[currPart] = value;
		else if (param == String("diffusionGridFile"))
			partDiffusionGridFile[currPart] = value;
		else if (param == String("diffusion"))
			part[currPart].diffusion = (float) strtod(value.val(), NULL);
		else if (param == String("charge"))
			part[currPart].charge = (float) strtod(value.val(), NULL);
		else if (param == String("radius"))
			part[currPart].radius = (float) strtod(value.val(), NULL);
		else if (param == String("eps"))
			part[currPart].eps = (float) strtod(value.val(), NULL);
		else if (param == String("reservoirFile"))
			partReservoirFile[currPart] = value;
		else if (param == String("tabulatedPotential"))
			tabulatedPotential = atoi(value.val());
		else if (param == String("tabulatedFile"))
			readTableFile(value, ++currTab);
		else if (param == String("tabulatedBondFile")) {
			if (numTabBondFiles >= btfcap) {
				String* temp = bondTableFile;
				btfcap *= 2;	
				bondTableFile = new String[btfcap];
				for (int j = 0; j < numTabBondFiles; j++)
					bondTableFile[j] = temp[j];
				delete[] temp;
			}
			if (readBondFile(value, ++currBond))
				numTabBondFiles++;
		} else if (param == String("inputParticles")) {
			if (readPartsFromFile) {
				printf("WARNING: More than one particle file specified. Ignoring new file.\n");
			} else {
				partFile = value;
				readPartsFromFile = true;
				loadedCoordinates = true;
			}
		} else if (param == String("inputBonds")) {
			if (readBondsFromFile) {
				printf("WARNING: More than one bond file specified. Ignoring new bond file.\n");
			} else {
				bondFile = value;				
				readBondsFromFile = true;
			}
		} else if (param == String("inputExcludes")) {
			if (readExcludesFromFile) {
				printf("WARNING: More than one exclude file specified. Ignoring new exclude file.\n");
			} else {
				excludeFile = value;				
				readExcludesFromFile = true;
			}
		} else if (param == String("exclude") or param == String("exclusion")) {
			excludeRule = value; 
		} else if (param == String("inputAngles")) {
			if (readAnglesFromFile) {
				printf("WARNING: More than one angle file specified. Ignoring new angle file.\n");
			} else {
				angleFile = value;
				readAnglesFromFile = true;
			}
		} else if (param == String("tabulatedAngleFile")) {
			if (numTabAngleFiles >= atfcap) {
				String* temp = angleTableFile;
				atfcap *= 2;	
				angleTableFile = new String[atfcap];
				for (int j = 0; j < numTabAngleFiles; j++)
					angleTableFile[j] = temp[j];
				delete[] temp;
			}
			if (readAngleFile(value, ++currAngle))
				numTabAngleFiles++;
		} else if (param == String("inputDihedrals")) {
			if (readDihedralsFromFile) {
				printf("WARNING: More than one dihedral file specified. Ignoring new dihedral file.\n");
			} else {
				dihedralFile = value;
				readDihedralsFromFile = true;
			}
		} else if (param == String("tabulatedDihedralFile")) {
			if (numTabDihedralFiles >= dtfcap) {
				String * temp = dihedralTableFile;
				dtfcap *= 2;
				dihedralTableFile = new String[dtfcap];
				for (int j = 0; j < numTabDihedralFiles; j++)
					dihedralTableFile[j] = temp[j];
				delete[] temp;
			}
			if (readDihedralFile(value, ++currDihedral))
				numTabDihedralFiles++;
		} else if (param == String("inputRestraints")) {
			if (readRestraintsFromFile) {
				printf("WARNING: More than one restraint file specified. Ignoring new restraint file.\n");
			} else {
				restraintFile = value;
				readRestraintsFromFile = true;
			}
		} else if (param == String("gridFileScale")) {
			partGridFileScale[currPart] = (float) strtod(value.val(), NULL);
		} else if (param == String("rigidBodyPotential")) {
			partRigidBodyGrid[currPart].push_back(value);
		}
                //Han-Yi Chou initial COM velocity for total particles
                else if (param == String("COM_Velocity"))
                    COM_Velocity = stringToVector3(value);

		// RIGID BODY
		else if (param == String("rigidBody")) {
			// part[++currPart] = BrownianParticleType(value);
			rigidBody[++currRB] = RigidBodyType(value, this);
			currPartClass = partClassRB;
		}
		else if (param == String("inertia"))
			rigidBody[currRB].inertia = stringToVector3( value );
		else if (param == String("rotDamping"))
			rigidBody[currRB].rotDamping = stringToVector3( value );

		else if (param == String("densityGrid"))
			rigidBody[currRB].addDensityGrid(value);
		else if (param == String("potentialGrid"))
			rigidBody[currRB].addPotentialGrid(value);
		else if (param == String("densityGridScale"))
			rigidBody[currRB].scaleDensityGrid(value);
		else if (param == String("potentialGridScale"))
			rigidBody[currRB].scalePotentialGrid(value);
		else if (param == String("pmfScale"))
			rigidBody[currRB].scalePMF(value);
		else if (param == String("position"))
			rigidBody[currRB].initPos = stringToVector3( value );
		else if (param == String("orientation"))
			rigidBody[currRB].initRot = stringToMatrix3( value );
                else if (param == String("momentum"))
                        rigidBody[currRB].initMomentum = stringToVector3(value);
                else if (param == String("angularMomentum"))//for angular momentum, serve as restart purpose
                        rigidBody[currRB].initAngularMomentum = stringToVector3(value);
		else if (param == String("inputRBCoordinates"))
			inputRBCoordinates = value;
		
		// COMMON
		else if (param == String("num")) {
			if (currPartClass == partClassPart)
				part[currPart].num = atoi(value.val());
			else if (currPartClass == partClassRB) 
				rigidBody[currRB].num = atoi(value.val());
		}
                //set mass here Han-Yi Chou
                else if (param == String("mass"))
                {
                    if(currPartClass == partClassPart)
                        part[currPart].mass    = (float) strtod(value.val(),NULL);
                    else if (currPartClass == partClassRB)
                        rigidBody[currRB].mass = (float) strtod(value.val(),NULL);
                }
                //set damping here, using anisotropic damping, i.e. data type Vector3 Han-Yi Chou
                else if (param == String("transDamping"))
                {
                    if(currPartClass == partClassPart)
                        part[currPart].transDamping    = stringToVector3(value);
                    else if (currPartClass == partClassRB)
                        rigidBody[currRB].transDamping = stringToVector3(value);

                }
		else if (param == String("gridFile")) {
			if (currPartClass == partClassPart)
				partGridFile[currPart] = value;
			else if (currPartClass == partClassRB)
				rigidBody[currRB].addPMF(value);
		}
		// UNKNOWN
		else {
			printf("ERROR: Unrecognized keyword `%s'.\n", param.val());
			exit(1);
		}
	}

	// extra configuration for RB types
	for (int i = 0; i < numRigidTypes; i++)
		rigidBody[i].setDampingCoeffs(timestep);