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 <iostream>
using namespace std;

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


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