Configuration.cpp

#include "Configuration.h"
#include "Angle.h"
#include "Dihedral.h"
#include "Restraint.h"
#include "ProductPotential.h"
#include <cmath>
#include <cassert>
#include <stdlib.h>     /* srand, rand */
#include <time.h>       /* time */
#include <string>
#include <iostream>
using namespace std;

#ifndef gpuErrchk
#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);
   }
}
#endif

namespace
{
    template<class T> 
    void convertString(const String& token, void* data)
    {
        exit(1);
    }

    template<> 
    void convertString<float>(const String& token, void* data)
    {
        float* tmp = (float*)data;
        *tmp = atof(token);
    }

    template<>
    void convertString<String>(const String& token, void* data)
    {
        String* tmp = (String*)data;
        *tmp = token;
    }

    template<class T>
    void stringToArray(String* str, int& size, T** array)
    {
        register int num;
        String *token;
        num =  str->tokenCount();
        size = num;
        *array = new T[num];
        token  = new String[num];
        str->tokenize(token);

        for(int i = 0; i < num; ++i)
            convertString<T>(token[i], (*array)+i);
        delete [] token;
    }
}
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());
		if (copyReplicaCoordinates <= 0) {
		    num /= simNum;
		}
  } 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

	// Count particles associated with rigid bodies
	num_rb_attached_particles = 0;
	if (numRigidTypes > 0) {
	    // grow list of rbs
	    for (int i = 0; i < numRigidTypes; i++) {
		RigidBodyType &rbt = rigidBody[i];
		rbt.attach_particles();
		num_rb_attached_particles += rbt.num * rbt.num_attached_particles();
	    }
	}
	assert( num_rb_attached_particles == 0 || simNum == 1 ); // replicas not yet implemented
	// num = num+num_rb_attached_particles;


	// Set the number capacity
	printf("\n%d particles\n", num);
	printf("%d particles attached to RBs\n", num_rb_attached_particles);

	if (numCap <= 0) numCap = numCapFactor*num; // max number of particles
	if (numCap <= 0) numCap = 20;

	if (readGroupSitesFromFile) readGroups();
	printf("%d groups\n", numGroupSites);

	// Allocate particle variables.
	// Each replica works with num+num_rb_attached_particles in array
	pos = new Vector3[ (num+num_rb_attached_particles) * simNum];

        //Han-Yi Chou
        if (ParticleDynamicType == String("Langevin") || ParticleDynamicType == String("NoseHooverLangevin"))
            momentum = new Vector3[(num+num_rb_attached_particles) * simNum];

	type   = new int[(num+num_rb_attached_particles) * simNum];
	serial = new int[(num+num_rb_attached_particles) * simNum];
	posLast = new Vector3[(num+num_rb_attached_particles) * simNum];

	{
	    int pidx = 0;
	    for (int i = 0; i < numRigidTypes; i++) { // Loop over RB types
		RigidBodyType &rbt = rigidBody[i];
		for (int j = 0; j < rbt.num; ++j) { // Loop over RBs
		    for (const int& t: rbt.get_attached_particle_types()) {
			type[num+pidx] = t;
			serial[num+pidx] = num+pidx;
			pidx++;
		    }
		}
	    }
	}	
	
        //Han-Yi Chou
        if(ParticleDynamicType == String("Langevin") || ParticleDynamicType == String("NoseHooverLangevin"))
           momLast = new Vector3[(num+num_rb_attached_particles) * simNum];
	name = new String[(num+num_rb_attached_particles) * 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 = find_particle_type(tokenList[2]);
				if (currType == -1) {
				    printf("Error: Unable to find particle type %s\n", tokenList[2].val());
				    exit(1);

				}
				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 (readBondAnglesFromFile) readBondAngles();
	if (readProductPotentialsFromFile) readProductPotentials();


	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);

	printf("Loading the potential grids...\n");
	// First load a single copy of each grid
	for (int i = 0; i < numParts; i++) 
        {
	    for(int j = 0; j < part[i].numPartGridFiles; ++j)
	    {
		std::string fname(partGridFile[i][j].val(), partGridFile[i][j].length());
		if (part_grid_dictionary.count( fname ) == 0)
		{
		    int len = fname.length();
		    if (len >= 3 && fname[len-3]=='.' && fname[len-2]=='d' && fname[len-1]=='x')
		    {
			part_grid_dictionary.insert({fname, BaseGrid(fname.c_str())});
		    }
		    else if  (len >= 4 && fname[len-4]=='.' && fname[len-3]=='d' && fname[len-2]=='e' && fname[len-1]=='f')
		    {
			assert(1==2); // Throw exception because this implementation needs to be revisited
/*                OverlordGrid* over = new OverlordGrid[part[i].numPartGridFiles];
		  part[i].meanPmf = new float[part[i].numPartGridFiles];
		  for(int j = 0; j < part[i].numPartGridFiles; ++j)
		  {
		  map = partGridFile[i][j];
		  len = map.length();
		  if (!(len >= 4 && map[len-4]=='.' && map[len-3]=='d' && map[len-2]=='e' && map[len-1]=='f'))
		  {
		  cout << "currently do not support different format " << endl;
		  exit(1);
		  }

		  String rootGrid = OverlordGrid::readDefFirst(map);
		  over[j] = 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].meanPmf[j] = part[i].pmf[j].mean();
		  }
		  part[i].pmf = static_cast<BaseGrid*>(over);
*/
		    } else {
			printf("WARNING: Unrecognized gridFile extension. Must be *.def or *.dx.\n");
			exit(-1);
		    }
		}
	    }
	}

	// Then assign grid addresses to particles
	for (int i = 0; i < numParts; i++)
        {
	    part[i].pmf     = new BaseGrid*[part[i].numPartGridFiles];
	    part[i].pmf_scale = new float[part[i].numPartGridFiles];
	    part[i].meanPmf = new float[part[i].numPartGridFiles];
	    for(int j = 0; j < part[i].numPartGridFiles; ++j)
	    {
		part[i].pmf[j] = &(part_grid_dictionary.find( std::string(partGridFile[i][j]) )->second);
		part[i].pmf_scale[j] = partGridFileScale[i][j];
		part[i].meanPmf[j] = part[i].pmf[j]->mean(); // TODO: review how this is used and decide whether to scale
	    }
		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);
	}

	{ // Add exclusions for RB attached particles
	    std::vector<Exclude> ex;
	    int start = num;
	    for (int i = 0; i < numRigidTypes; i++) { // Loop over RB types
		RigidBodyType &rbt = rigidBody[i];
		const int nap = rbt.num_attached_particles();
		for (int j = 0; j < rbt.num; ++j) { // Loop over RBs
		    for (int ai = 0; ai < nap-1; ++ai) {
			for (int aj = ai+1; aj < nap; ++aj) {
			    ex.push_back( Exclude( ai+start, aj+start ) );
			}
		    }
		    start += nap;
		}
	    }
	    // copy
	    int oldNumExcludes = numExcludes;
	    numExcludes = numExcludes + ex.size();
	    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] = ex[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+num_rb_attached_particles; ++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;
	for(int i = 0; i < numParts; ++i)
        {
            if(partGridFile[i] != NULL) 
            {
                delete[] partGridFile[i];
                partGridFile[i] = NULL;
            }
            if(partGridFileScale[i] != NULL)
            {
                delete[] partGridFileScale[i];
                partGridFileScale[i] = NULL;
            }
        }
        delete [] partGridFile;
        delete [] partGridFileScale;
        //delete numPartGridFiles;
	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;
	if (bondAngles != NULL) delete[] bondAngles;
	if (productPotentials != NULL) delete[] productPotentials;

	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 grids to GPU %d\n", GPUManager::current());
    for (const auto& pair : part_grid_dictionary)
    {
	// Copy PMF
	const BaseGrid& g = pair.second;
	BaseGrid *g_d = g.copy_to_cuda();
	part_grid_dictionary_d.insert({pair.first, g_d});
	// printf("Assigning grid for %s to %p (originally %p)\n", pair.first.c_str(), (void *) part_grid_dictionary_d[pair.first], (void *) g_d);
    }

    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++) 
        {
		BrownianParticleType *b = new BrownianParticleType(part[i]);
		// Copy PMF pointers
		if (part[i].pmf != NULL) 
                {
		    {
			BaseGrid** tmp_d = new BaseGrid*[part[i].numPartGridFiles];
			BaseGrid** tmp   = new BaseGrid*[part[i].numPartGridFiles];
			for(int j = 0; j < part[i].numPartGridFiles; ++j) {
			    // printf("Retrieving grid for %s (at %p)\n", partGridFile[i][j].val(), (void *) part_grid_dictionary_d[std::string(partGridFile[i][j])]);
			    tmp[j] = part_grid_dictionary_d[std::string(partGridFile[i][j])];
			}
			gpuErrchk(cudaMalloc(&tmp_d, sizeof(BaseGrid*)*part[i].numPartGridFiles));
			gpuErrchk(cudaMemcpy(tmp_d, tmp, sizeof(BaseGrid*)*part[i].numPartGridFiles,
					     cudaMemcpyHostToDevice));
			b->pmf = tmp_d;
		    }

		    {
			float *tmp;
			gpuErrchk(cudaMalloc(&tmp, sizeof(float)*part[i].numPartGridFiles));
			gpuErrchk(cudaMemcpy(tmp, part[i].pmf_scale, sizeof(float)*part[i].numPartGridFiles,
					     cudaMemcpyHostToDevice));
			b->pmf_scale = tmp;
		    }

		    {
			float *tmp;
			gpuErrchk(cudaMalloc(&tmp, sizeof(float)*part[i].numPartGridFiles));
			gpuErrchk(cudaMemcpy(tmp, part[i].meanPmf, sizeof(float)*part[i].numPartGridFiles, 
					     cudaMemcpyHostToDevice));
			b->meanPmf = tmp;
		    }

		    {
			BoundaryCondition *tmp;
			size_t s = sizeof(BoundaryCondition)*part[i].numPartGridFiles;
			gpuErrchk(cudaMalloc(&tmp, s));
			gpuErrchk(cudaMemcpy(tmp, part[i].pmf_boundary_conditions, s, cudaMemcpyHostToDevice));
			b->pmf_boundary_conditions = tmp;
		    }
                    
		}

		// Copy the diffusion grid
		if (part[i].diffusionGrid != NULL) {
		    b->diffusionGrid = part[i].diffusionGrid->copy_to_cuda();
		} else {
		    b->diffusionGrid = NULL;
		}
		
		//b->pmf = pmf;
		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_rb_attached_particles) * 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;

	unsigned long int r0 = clock();
	for (int i = 0; i < 4; i++)
	    r0 *= r0 + 1;
	seed = time(NULL) + r0;

	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 = 0;
	//	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;
	numGroupSites = 0;
	readGroupSitesFromFile = 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;

	readBondAnglesFromFile = false;
	numBondAngles = 0;
	bondAngles = NULL;

	readProductPotentialsFromFile = false;
	numProductPotentials = 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;

        ParticleInterpolationType = 0;
        RigidBodyInterpolationType = 0;
}

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];
	partGridFile       = new String*[numParts];
        //partGridFileScale = new float[numParts];
        partGridFileScale  = new float*[numParts];
        //int numPartGridFiles = new int[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;
	}*/

        for(int i = 0; i < numParts; ++i)
        {
            partGridFile[i] = NULL;
            partGridFileScale[i] = NULL;
            //part[i].numPartGridFiles = -1;
        }
        //for(int i = 0; i < numParts; ++i)
          //  cout << part[i].numPartGridFiles << endl;

	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);
		// printf("Parsing %s: %s\n", param.val(), value.val());
		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());