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

	{
	    for (size_t replica = 0; replica < simNum; ++replica) {
		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()) {
			    size_t _idx = num+pidx+replica*(num+num_rb_attached_particles);
			    type[_idx] = t;
			    serial[_idx] = 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");

		// Copy restart data to replicas
		for (size_t replica = 1; replica < simNum; ++replica) {
		    for (size_t i = 0; i < num; ++i) {
			size_t _idx = i+replica*(num+num_rb_attached_particles);
			type[_idx] = type[i];
			serial[_idx] = serial[i];
			pos[_idx] = pos[i];
		    }
		}

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