diff --git a/src/ComputeForce.cu b/src/ComputeForce.cu index 4efa1dbb24b7b1ea674523a0f98a6b2fd3019031..a7c8dd8bc5989a77d0aa2d76f3a48c1c82586ddc 100644 --- a/src/ComputeForce.cu +++ b/src/ComputeForce.cu @@ -43,7 +43,7 @@ ComputeForce::ComputeForce(const Configuration& c, const int numReplicas = 1) : numExcludes(c.numExcludes), numAngles(c.numAngles), numTabAngleFiles(c.numTabAngleFiles), numDihedrals(c.numDihedrals), numTabDihedralFiles(c.numTabDihedralFiles), numRestraints(c.numRestraints), - numBondAngles(c.numBondAngles), numCrossPotentials(c.numCrossPotentials), + numBondAngles(c.numBondAngles), numProductPotentials(c.numProductPotentials), numReplicas(numReplicas) { // Allocate the parameter tables. @@ -140,7 +140,7 @@ ComputeForce::ComputeForce(const Configuration& c, const int numReplicas = 1) : restraintIds_d = NULL; bondAngleList_d = NULL; - cross_potential_list_d = NULL; + product_potential_list_d = NULL; //Calculate the number of blocks the grid should contain gridSize = num / NUM_THREADS + 1; @@ -746,9 +746,9 @@ float ComputeForce::computeTabulated(bool get_energy) { //Mlog: the commented function doesn't use bondList, uncomment for testing. //if(bondMap_d != NULL && tableBond_d != NULL) - if(cross_potential_list_d != NULL && cross_potentials_d != NULL) + if(product_potential_list_d != NULL && product_potentials_d != NULL) { - computeCrossPotentials <<<nb, numThreads, 0, gpuman.get_next_stream()>>> ( forceInternal_d, pos_d, sys_d, numReplicas*numCrossPotentials, cross_potential_particles_d, cross_potentials_d, cross_potential_list_d, numCrossed_d, energies_d, get_energy); + computeProductPotentials <<<nb, numThreads, 0, gpuman.get_next_stream()>>> ( forceInternal_d, pos_d, sys_d, numReplicas*numProductPotentials, product_potential_particles_d, product_potentials_d, product_potential_list_d, productCount_d, energies_d, get_energy); } if(bondAngleList_d != NULL && tableBond_d != NULL && tableAngle_d != NULL) @@ -903,7 +903,7 @@ void ComputeForce::setForceInternalOnDevice(Vector3* f) { } -void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap, Angle* angles, Dihedral* dihedrals, const Restraint* const restraints, const BondAngle* const bondAngles, const XpotMap simple_potential_map, const std::vector<SimplePotential> simple_potentials, const CrossPotentialConf* const cross_potential_confs) +void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap, Angle* angles, Dihedral* dihedrals, const Restraint* const restraints, const BondAngle* const bondAngles, const XpotMap simple_potential_map, const std::vector<SimplePotential> simple_potentials, const ProductPotentialConf* const product_potential_confs) { // type_d gpuErrchk(cudaMalloc(&type_d, sizeof(int) * num * simNum)); @@ -1002,66 +1002,66 @@ void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap, } - if (numCrossPotentials > 0) { + if (numProductPotentials > 0) { // Count particles int n_pots = 0; int n_particles = 0; - for (int i=0; i < numCrossPotentials; ++i) { - const CrossPotentialConf& c = cross_potential_confs[i]; + for (int i=0; i < numProductPotentials; ++i) { + const ProductPotentialConf& c = product_potential_confs[i]; n_pots += c.indices.size(); for (int j=0; j < c.indices.size(); ++j) { n_particles += c.indices[j].size(); } } - // printf("DEBUG: Found %d particles participating in %d potentials forming %d crossPotentials\n", - // n_particles, n_pots, numCrossPotentials); + // printf("DEBUG: Found %d particles participating in %d potentials forming %d productPotentials\n", + // n_particles, n_pots, numProductPotentials); - // Build crossPotentialLists on host + // Build productPotentialLists on host int *particle_list = new int[n_particles]; - SimplePotential *cross_potentials = new SimplePotential[n_pots]; - uint2 *cross_potential_list = new uint2[numCrossPotentials]; - unsigned short *numCrossed = new unsigned short[numCrossPotentials]; + SimplePotential *product_potentials = new SimplePotential[n_pots]; + uint2 *product_potential_list = new uint2[numProductPotentials]; + unsigned short *productCount = new unsigned short[numProductPotentials]; n_particles = 0; n_pots = 0; - for (int i=0; i < numCrossPotentials; ++i) { - const CrossPotentialConf& c = cross_potential_confs[i]; - cross_potential_list[i] = make_uint2( n_pots, n_particles ); + for (int i=0; i < numProductPotentials; ++i) { + const ProductPotentialConf& c = product_potential_confs[i]; + product_potential_list[i] = make_uint2( n_pots, n_particles ); for (int j=0; j < c.indices.size(); ++j) { unsigned int sp_i = simple_potential_map.at(c.potential_names[j]); - cross_potentials[n_pots] = simple_potentials[sp_i]; - cross_potentials[n_pots++].pot = simple_potential_pots_d[sp_i]; + product_potentials[n_pots] = simple_potentials[sp_i]; + product_potentials[n_pots++].pot = simple_potential_pots_d[sp_i]; for (int k=0; k < c.indices[j].size(); ++k) { particle_list[n_particles++] = c.indices[j][k]; } } - numCrossed[i] = c.indices.size(); + productCount[i] = c.indices.size(); } // Copy to device size_t sz = sizeof(int)*n_particles; - gpuErrchk(cudaMalloc(&cross_potential_particles_d, sz)); - gpuErrchk(cudaMemcpyAsync(cross_potential_particles_d, particle_list, sz, + gpuErrchk(cudaMalloc(&product_potential_particles_d, sz)); + gpuErrchk(cudaMemcpyAsync(product_potential_particles_d, particle_list, sz, cudaMemcpyHostToDevice)); sz = sizeof(SimplePotential)*n_pots; - gpuErrchk(cudaMalloc(&cross_potentials_d, sz)); - gpuErrchk(cudaMemcpyAsync(cross_potentials_d, cross_potentials, sz, + gpuErrchk(cudaMalloc(&product_potentials_d, sz)); + gpuErrchk(cudaMemcpyAsync(product_potentials_d, product_potentials, sz, cudaMemcpyHostToDevice)); - sz = sizeof(uint2)*numCrossPotentials; - gpuErrchk(cudaMalloc(&cross_potential_list_d, sz)); - gpuErrchk(cudaMemcpyAsync(cross_potential_list_d, cross_potential_list, sz, + sz = sizeof(uint2)*numProductPotentials; + gpuErrchk(cudaMalloc(&product_potential_list_d, sz)); + gpuErrchk(cudaMemcpyAsync(product_potential_list_d, product_potential_list, sz, cudaMemcpyHostToDevice)); - sz = sizeof(unsigned short)*numCrossPotentials; - gpuErrchk(cudaMalloc(&numCrossed_d, sz)); - gpuErrchk(cudaMemcpyAsync(numCrossed_d, numCrossed, sz, + sz = sizeof(unsigned short)*numProductPotentials; + gpuErrchk(cudaMalloc(&productCount_d, sz)); + gpuErrchk(cudaMemcpyAsync(productCount_d, productCount, sz, cudaMemcpyHostToDevice)); // Clean up delete[] particle_list; - delete[] cross_potentials; - delete[] cross_potential_list; - delete[] numCrossed; + delete[] product_potentials; + delete[] product_potential_list; + delete[] productCount; } gpuErrchk(cudaDeviceSynchronize()); diff --git a/src/ComputeForce.cuh b/src/ComputeForce.cuh index 07c023a7ed6fab223773d8b4122df33440fe109e..268ec05c4f7902d0f59145c5bb93db07b64bc8fb 100644 --- a/src/ComputeForce.cuh +++ b/src/ComputeForce.cuh @@ -1094,24 +1094,26 @@ void computeTabulatedBondAngles(Vector3* force, } __global__ -void computeCrossPotentials(Vector3* force, - Vector3* __restrict__ pos, - BaseGrid* __restrict__ sys, - int numCrossPotentials, - int* __restrict__ crossPotentialParticles, - SimplePotential* __restrict__ potentialList, - uint2* __restrict__ crossPotential_list, - unsigned short* __restrict__ numCrossed, - float* energy, bool get_energy) { +void computeProductPotentials(Vector3* force, + Vector3* __restrict__ pos, + BaseGrid* __restrict__ sys, + int numProductPotentials, + int* __restrict__ productPotentialParticles, + SimplePotential* __restrict__ potentialList, + uint2* __restrict__ productPotential_list, + unsigned short* __restrict__ productCount, + float* energy, bool get_energy) { /* - crossPotential_list[i].x : index of first potential in potentialList for i_th crossPotential - crossPotential_list[i].y : index of first atom in crossPotentialParticles for i_th crossPotential + productPotential_list[i].x : index of first potential in potentialList for i_th productPotential + productPotential_list[i].y : index of first atom in productPotentialParticles for i_th productPotential - for three potentials, angle, bond, angle, we would have the following atomic indices in crossPotentialParticles: - pot1 : crossPotential_list[i].y, crossPotential_list[i].y + 1 , crossPotential_list[i].y + 2 - pot2 : crossPotential_list[i].y + 3, crossPotential_list[i].y + 4 + for three potentials, angle, bond, angle, we would have the following atomic indices in productPotentialParticles: + pot1 : productPotential_list[i].y, productPotential_list[i].y + 1 , productPotential_list[i].y + 2 + pot2 : productPotential_list[i].y + 3, productPotential_list[i].y + 4 and - pot3 : crossPotential_list[i].y + 5, crossPotential_list[i].y + 6, crossPotential_list[i].y + 7 + pot3 : productPotential_list[i].y + 5, productPotential_list[i].y + 6, productPotential_list[i].y + 7 + + productCount[i] : number of potentials in the i_th productPotential */ // CRAPPY NAIVE IMPLEMENTATION @@ -1119,22 +1121,22 @@ void computeCrossPotentials(Vector3* force, float2 energy_and_deriv[MAX_XPOTS]; float tmp_force; - for (int i = threadIdx.x+blockIdx.x*blockDim.x; i<numCrossPotentials; i+=blockDim.x*gridDim.x) { - unsigned short num_pots = numCrossed[i]; + for (int i = threadIdx.x+blockIdx.x*blockDim.x; i<numProductPotentials; i+=blockDim.x*gridDim.x) { + unsigned short num_pots = productCount[i]; - unsigned int part_idx = crossPotential_list[i].y; + unsigned int part_idx = productPotential_list[i].y; #pragma unroll for (unsigned short int j = 0; j < MAX_XPOTS; ++j) { if (j == num_pots) break; - SimplePotential& p = potentialList[ crossPotential_list[i].x + j ]; + SimplePotential& p = potentialList[ productPotential_list[i].x + j ]; // Hidden branch divergence in compute_value => sort potentials by type before running kernel - float tmp = p.compute_value(pos,sys, &crossPotentialParticles[part_idx]); + float tmp = p.compute_value(pos,sys, &productPotentialParticles[part_idx]); energy_and_deriv[j] = p.compute_energy_and_deriv(tmp); part_idx += p.type==BOND? 2: p.type==ANGLE? 3: 4; } - part_idx = crossPotential_list[i].y; + part_idx = productPotential_list[i].y; #pragma unroll for (unsigned short int j = 0; j < MAX_XPOTS; ++j) { if (j == num_pots) break; @@ -1149,8 +1151,8 @@ void computeCrossPotentials(Vector3* force, // TODO add energy // TODO make it work with replicas - SimplePotential& p = potentialList[ crossPotential_list[i].x + j ]; - p.apply_force(pos,sys, force, &crossPotentialParticles[part_idx], tmp_force); + SimplePotential& p = potentialList[ productPotential_list[i].x + j ]; + p.apply_force(pos,sys, force, &productPotentialParticles[part_idx], tmp_force); part_idx += p.type==BOND? 2: p.type==ANGLE? 3: 4; } } diff --git a/src/ComputeForce.h b/src/ComputeForce.h index a1e142638d08f34768ff9541e7676ecc01e38dde..8270fd5f6062ade18bc4ba80942ea12576ff2f01 100644 --- a/src/ComputeForce.h +++ b/src/ComputeForce.h @@ -24,7 +24,7 @@ #include "TabulatedPotential.h" #include "TabulatedAngle.h" #include "TabulatedDihedral.h" -#include "CrossPotential.h" +#include "ProductPotential.h" #include "GPUManager.h" // #include <map> @@ -90,7 +90,7 @@ public: void copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap, Exclude* excludes, int2* excludeMap, Angle* angles, Dihedral* dihedrals, const Restraint* const restraints, const BondAngle* const bondAngles, const XpotMap simple_potential_map, const std::vector<SimplePotential> simple_potentials, - const CrossPotentialConf* const cross_potential_confs); + const ProductPotentialConf* const product_potential_confs); void copyToCUDA(Vector3* forceInternal, Vector3* pos, Vector3* mom); void copyToCUDA(Vector3* forceInternal, Vector3* pos, Vector3* mom, float* random); @@ -239,14 +239,13 @@ private: BondAngle* bondAngles_d; int4* bondAngleList_d; - int numCrossPotentials; + int numProductPotentials; float** simple_potential_pots_d; SimplePotential* simple_potentials_d; - int* cross_potential_particles_d; - SimplePotential* cross_potentials_d; - uint2* cross_potential_list_d; - // ushort4* cross_potential_particle_counts_d; - unsigned short* numCrossed_d; + int* product_potential_particles_d; + SimplePotential* product_potentials_d; + uint2* product_potential_list_d; + unsigned short* productCount_d; int3* bondList_d; int4* angleList_d; diff --git a/src/Configuration.cpp b/src/Configuration.cpp index ddd26fbda8099ceb67b8853986b59f48a7ca4d9b..89ac7c0333e60527e93bfd97a2eadc03ed8d97de 100644 --- a/src/Configuration.cpp +++ b/src/Configuration.cpp @@ -2,7 +2,7 @@ #include "Angle.h" #include "Dihedral.h" #include "Restraint.h" -#include "CrossPotential.h" +#include "ProductPotential.h" #include <cmath> #include <cassert> #include <stdlib.h> /* srand, rand */ @@ -229,7 +229,7 @@ Configuration::Configuration(const char* config_file, int simNum, bool debug) : if (readDihedralsFromFile) readDihedrals(); if (readRestraintsFromFile) readRestraints(); if (readBondAnglesFromFile) readBondAngles(); - if (readCrossPotentialsFromFile) readCrossPotentials(); + if (readProductPotentialsFromFile) readProductPotentials(); if (temperatureGridFile.length() != 0) { @@ -565,7 +565,7 @@ Configuration::~Configuration() { if (angles != NULL) delete[] angles; if (dihedrals != NULL) delete[] dihedrals; if (bondAngles != NULL) delete[] bondAngles; - if (crossPotentials != NULL) delete[] crossPotentials; + if (productPotentials != NULL) delete[] productPotentials; delete[] numPartsOfType; @@ -790,8 +790,8 @@ void Configuration::setDefaults() { numBondAngles = 0; bondAngles = NULL; - readCrossPotentialsFromFile = false; - numCrossPotentials = 0; + readProductPotentialsFromFile = false; + numProductPotentials = 0; readRestraintsFromFile = false; @@ -1041,12 +1041,12 @@ int Configuration::readParameters(const char * config_file) { bondAngleFile = value; readBondAnglesFromFile = true; } - } else if (param == String("inputCrossPotentials")) { + } else if (param == String("inputProductPotentials")) { if (readBondAnglesFromFile) { - printf("WARNING: More than one cross potential file specified. Ignoring new file.\n"); + printf("WARNING: More than one product potential file specified. Ignoring new file.\n"); } else { - crossPotentialFile = value; - readCrossPotentialsFromFile = true; + productPotentialFile = value; + readProductPotentialsFromFile = true; } } else if (param == String("tabulatedAngleFile")) { if (numTabAngleFiles >= atfcap) { @@ -1823,20 +1823,20 @@ void Configuration::readBondAngles() { // angles[i].print(); } -void Configuration::readCrossPotentials() { - FILE* inp = fopen(crossPotentialFile.val(), "r"); +void Configuration::readProductPotentials() { + FILE* inp = fopen(productPotentialFile.val(), "r"); char line[256]; int capacity = 256; - numCrossPotentials = 0; - crossPotentials = new CrossPotentialConf[capacity]; + numProductPotentials = 0; + productPotentials = new ProductPotentialConf[capacity]; // If the angle file cannot be found, exit the program if (inp == NULL) { - printf("WARNING: Could not open `%s'.\n", crossPotentialFile.val()); - printf("This simulation will not use cross potentials.\n"); + printf("WARNING: Could not open `%s'.\n", productPotentialFile.val()); + printf("This simulation will not use product potentials.\n"); return; } - printf("DEBUG: READING CROSSPOT FILE\n"); + printf("DEBUG: READING PRODUCT POTENTAL FILE\n"); std::vector<std::vector<int>> indices; std::vector<int> tmp; std::vector<String> pot_names; @@ -1852,12 +1852,12 @@ void Configuration::readCrossPotentials() { tmp.clear(); pot_names.clear(); - printf("\rDEBUG: reading line %d",numCrossPotentials+1); + printf("\rDEBUG: reading line %d",numProductPotentials+1); - // Legitimate CrossPotential inputs have at least 7 tokens + // Legitimate ProductPotential inputs have at least 7 tokens // BONDANGLE | INDEX1 | INDEX2 | INDEX3 | POT_FILENAME1 | INDEX4 | INDEX5 | POT_FILENAME2 ... if (numTokens < 7) { - printf("WARNING: Invalid cross potential input line (too few tokens %d): %s\n", numTokens, line); + printf("WARNING: Invalid product potential input line (too few tokens %d): %s\n", numTokens, line); continue; } @@ -1877,7 +1877,7 @@ void Configuration::readCrossPotentials() { // Could not find fileName in dictionary, so read and add it unsigned int s = tmp.size(); if (s < 2 || s > 4) { - printf("WARNING: Invalid cross potential input line (indices of potential %d == %d): %s\n", i, s, line); + printf("WARNING: Invalid product potential input line (indices of potential %d == %d): %s\n", i, s, line); continue; } @@ -1894,21 +1894,21 @@ void Configuration::readCrossPotentials() { } } - if (numCrossPotentials >= capacity) { - CrossPotentialConf* temp = crossPotentials; + if (numProductPotentials >= capacity) { + ProductPotentialConf* temp = productPotentials; capacity *= 2; - crossPotentials = new CrossPotentialConf[capacity]; - for (int i = 0; i < numCrossPotentials; i++) - crossPotentials[i] = temp[i]; + productPotentials = new ProductPotentialConf[capacity]; + for (int i = 0; i < numProductPotentials; i++) + productPotentials[i] = temp[i]; delete[] temp; } - CrossPotentialConf a(indices, pot_names); - crossPotentials[numCrossPotentials++] = a; + ProductPotentialConf a(indices, pot_names); + productPotentials[numProductPotentials++] = a; delete[] tokenList; } printf("\nDEBUG: Sorting\n"); - std::sort(crossPotentials, crossPotentials + numCrossPotentials, compare()); + std::sort(productPotentials, productPotentials + numProductPotentials, compare()); // for(int i = 0; i < numAngles; i++) // angles[i].print(); @@ -2528,7 +2528,7 @@ bool Configuration::compare::operator()(const BondAngle& lhs, const BondAngle& r return lhs.ind4 < rhs.ind4; } -bool Configuration::compare::operator()(const CrossPotentialConf& lhs, const CrossPotentialConf& rhs) { +bool Configuration::compare::operator()(const ProductPotentialConf& lhs, const ProductPotentialConf& rhs) { int diff = rhs.indices.size() - lhs.indices.size(); if (diff != 0) return diff > 0; diff --git a/src/Configuration.h b/src/Configuration.h index 64b7bb5dd66574270806367f51deb2c88ee61247..55b0ad2f11027a2a0e6cb6f5343ef34dd9b7f313 100644 --- a/src/Configuration.h +++ b/src/Configuration.h @@ -24,7 +24,7 @@ #include "TrajectoryWriter.h" #include "TabulatedPotential.h" #include "TabulatedAngle.h" -#include "CrossPotential.h" +#include "ProductPotential.h" #include "GPUManager.h" #include "RigidBodyType.h" #include "RigidBody.h" @@ -49,7 +49,7 @@ class Configuration { bool operator()(const Angle& lhs, const Angle& rhs); bool operator()(const Dihedral& lhs, const Dihedral& rhs); bool operator()(const BondAngle& lhs, const BondAngle& rhs); - bool operator()(const CrossPotentialConf& lhs, const CrossPotentialConf& rhs); + bool operator()(const ProductPotentialConf& lhs, const ProductPotentialConf& rhs); }; void setDefaults(); @@ -246,11 +246,11 @@ public: Restraint* restraints; - void readCrossPotentials(); - String crossPotentialFile; - int numCrossPotentials; - bool readCrossPotentialsFromFile; - CrossPotentialConf* crossPotentials; + void readProductPotentials(); + String productPotentialFile; + int numProductPotentials; + bool readProductPotentialsFromFile; + ProductPotentialConf* productPotentials; // boost::unordered_map<String, unsigned int> simple_potential_ids; XpotMap simple_potential_ids; std::vector<SimplePotential> simple_potentials; diff --git a/src/GrandBrownTown.cu b/src/GrandBrownTown.cu index 7b19f7c60145fa78ef2e8f9db4ea9a258a0973af..1c2a5630ec764347091c92f7aac1d8ce0280fa3e 100644 --- a/src/GrandBrownTown.cu +++ b/src/GrandBrownTown.cu @@ -254,7 +254,7 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg, internal = new ComputeForce(c, numReplicas); //MLog: I did the other halve of the copyToCUDA function from the Configuration class here, keep an eye on any mistakes that may occur due to the location. - internal -> copyToCUDA(c.simNum, c.type, c.bonds, c.bondMap, c.excludes, c.excludeMap, c.angles, c.dihedrals, c.restraints, c.bondAngles, c.simple_potential_ids, c.simple_potentials, c.crossPotentials ); + internal -> copyToCUDA(c.simNum, c.type, c.bonds, c.bondMap, c.excludes, c.excludeMap, c.angles, c.dihedrals, c.restraints, c.bondAngles, c.simple_potential_ids, c.simple_potentials, c.productPotentials ); // TODO: check for duplicate potentials if (c.tabulatedPotential) { diff --git a/src/ProductPotential.h b/src/ProductPotential.h new file mode 100644 index 0000000000000000000000000000000000000000..ff918e609c1b59827b71bfa0af68401ba287079a --- /dev/null +++ b/src/ProductPotential.h @@ -0,0 +1,30 @@ +#pragma once + +#ifdef __CUDACC__ + #define HOST __host__ + #define DEVICE __device__ +#else + #define HOST + #define DEVICE +#endif + +#include "useful.h" +#include <cuda.h> + + +class ProductPotentialConf { +public: + ProductPotentialConf() {} + ProductPotentialConf( std::vector< std::vector<int> > indices, std::vector<String> potential_names ) : + indices(indices), potential_names(potential_names) { } + + std::vector< std::vector<int> > indices; /* indices of particles */ + std::vector<String> potential_names; + + inline ProductPotentialConf(const ProductPotentialConf& a) : indices(a.indices), potential_names(a.potential_names) { } + + + /* String toString(); */ + /* void print(); */ +}; +