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