From 4e126a49b63b870bb619496157f20a1bbec40b48 Mon Sep 17 00:00:00 2001 From: Chris Maffeo <cmaffeo2@illinois.edu> Date: Thu, 6 Aug 2020 18:01:43 -0500 Subject: [PATCH] Changed name to product potential --- src/ComputeForce.cu | 67 +++++++++++++++++++++--------------------- src/ComputeForce.cuh | 48 +++++++++++++++--------------- src/ComputeForce.h | 15 +++++----- src/Configuration.cpp | 58 ++++++++++++++++++------------------ src/Configuration.h | 14 ++++----- src/GrandBrownTown.cu | 2 +- src/ProductPotential.h | 30 +++++++++++++++++++ 7 files changed, 132 insertions(+), 102 deletions(-) create mode 100644 src/ProductPotential.h diff --git a/src/ComputeForce.cu b/src/ComputeForce.cu index f5b42bc..253ff7e 100644 --- a/src/ComputeForce.cu +++ b/src/ComputeForce.cu @@ -49,8 +49,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), numGroupSites(c.numGroupSites), numReplicas(numReplicas) { @@ -236,7 +235,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_rb_attached_particles) / NUM_THREADS + 1; @@ -731,9 +730,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) @@ -892,7 +891,7 @@ void ComputeForce::setForceInternalOnDevice(Vector3* f) { gpuErrchk(cudaMemcpy(forceInternal_d[0], f, sizeof(Vector3) * tot_num, cudaMemcpyHostToDevice)); } -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) { assert(simNum == numReplicas); // Not sure why we have both of these things int tot_num_with_rb = (num+num_rb_attached_particles) * simNum; @@ -994,66 +993,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 eb01f2f..275764f 100644 --- a/src/ComputeForce.cuh +++ b/src/ComputeForce.cuh @@ -890,24 +890,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 @@ -915,22 +917,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; @@ -945,8 +947,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 1801267..d774e6a 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); @@ -268,14 +268,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 2ce832f..b08fdd3 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 */ @@ -273,7 +273,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) { @@ -635,7 +635,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; @@ -871,8 +871,8 @@ void Configuration::setDefaults() { numBondAngles = 0; bondAngles = NULL; - readCrossPotentialsFromFile = false; - numCrossPotentials = 0; + readProductPotentialsFromFile = false; + numProductPotentials = 0; readRestraintsFromFile = false; @@ -1143,12 +1143,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) { @@ -2016,20 +2016,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; @@ -2045,12 +2045,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; } @@ -2070,7 +2070,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; } @@ -2087,21 +2087,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(); @@ -2721,7 +2721,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 00fddec..f956b92 100644 --- a/src/Configuration.h +++ b/src/Configuration.h @@ -25,7 +25,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" @@ -50,7 +50,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(); @@ -267,11 +267,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 27e02b0..ef6ef5d 100644 --- a/src/GrandBrownTown.cu +++ b/src/GrandBrownTown.cu @@ -262,7 +262,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 ); if (numGroupSites > 0) init_cuda_group_sites(); // TODO: check for duplicate potentials diff --git a/src/ProductPotential.h b/src/ProductPotential.h new file mode 100644 index 0000000..ff918e6 --- /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(); */ +}; + -- GitLab