From 4c59c134345e74ad7c2b0b5d88d5402f724fa5ac Mon Sep 17 00:00:00 2001
From: Chris Maffeo <cmaffeo2@illinois.edu>
Date: Wed, 30 Dec 2020 13:38:03 -0600
Subject: [PATCH] Added support for replicas to productPotential

---
 src/ComputeForce.cu  | 44 +++++++++++++++++++++++++-------------------
 src/ComputeForce.cuh |  1 -
 2 files changed, 25 insertions(+), 20 deletions(-)

diff --git a/src/ComputeForce.cu b/src/ComputeForce.cu
index 253ff7e..af9a276 100644
--- a/src/ComputeForce.cu
+++ b/src/ComputeForce.cu
@@ -1008,42 +1008,48 @@ void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap,
 	    // 	   n_particles, n_pots, numProductPotentials);
 
 	    // Build productPotentialLists on host
-	    int *particle_list = new int[n_particles];
+	    int *particle_list = new int[n_particles*numReplicas];
 	    SimplePotential *product_potentials = new SimplePotential[n_pots];
-	    uint2 *product_potential_list = new uint2[numProductPotentials];
-	    unsigned short *productCount = new unsigned short[numProductPotentials];
+	    uint2 *product_potential_list = new uint2[numProductPotentials*numReplicas];
+	    unsigned short *productCount = new unsigned short[numProductPotentials*numReplicas];
 
 	    n_particles = 0;
-	    n_pots = 0;
-	    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]);
-		    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];
+	    
+	    for (unsigned int r=0; r < numReplicas; ++r) {
+		n_pots = 0;
+		for (int i=0; i < numProductPotentials; ++i) {
+		    const ProductPotentialConf& c = product_potential_confs[i];
+		    product_potential_list[i+r*numProductPotentials] = make_uint2( n_pots, n_particles );
+
+		    for (int j=0; j < c.indices.size(); ++j) {
+			if (r == 0) {
+			    unsigned int sp_i = simple_potential_map.at(c.potential_names[j]);
+			    product_potentials[n_pots] = simple_potentials[sp_i];
+			    product_potentials[n_pots].pot = simple_potential_pots_d[sp_i];
+			}
+			++n_pots;
+			for (int k=0; k < c.indices[j].size(); ++k) {
+			    particle_list[n_particles++] = c.indices[j][k]+r*num;
+			}
 		    }
+		    productCount[i+r*numProductPotentials] = c.indices.size();
 		}
-		productCount[i] = c.indices.size();
 	    }
 
 	    // Copy to device
-	    size_t sz = sizeof(int)*n_particles;
+	    size_t sz = n_particles*numReplicas * sizeof(int);
 	    gpuErrchk(cudaMalloc(&product_potential_particles_d, sz));
 	    gpuErrchk(cudaMemcpyAsync(product_potential_particles_d, particle_list, sz,
 	    				  cudaMemcpyHostToDevice));
-	    sz = sizeof(SimplePotential)*n_pots;
+	    sz = n_pots * sizeof(SimplePotential);
 	    gpuErrchk(cudaMalloc(&product_potentials_d, sz));
 	    gpuErrchk(cudaMemcpyAsync(product_potentials_d, product_potentials, sz,
 	    				  cudaMemcpyHostToDevice));
-	    sz = sizeof(uint2)*numProductPotentials;
+	    sz = numProductPotentials*numReplicas * sizeof(uint2);
 	    gpuErrchk(cudaMalloc(&product_potential_list_d, sz));
 	    gpuErrchk(cudaMemcpyAsync(product_potential_list_d, product_potential_list, sz,
 	    				  cudaMemcpyHostToDevice));
-	    sz = sizeof(unsigned short)*numProductPotentials;
+	    sz = numProductPotentials*numReplicas * sizeof(unsigned short);
 	    gpuErrchk(cudaMalloc(&productCount_d, sz));
 	    gpuErrchk(cudaMemcpyAsync(productCount_d, productCount, sz,
 	    				  cudaMemcpyHostToDevice));
diff --git a/src/ComputeForce.cuh b/src/ComputeForce.cuh
index c14c8b7..5d9cfb4 100644
--- a/src/ComputeForce.cuh
+++ b/src/ComputeForce.cuh
@@ -946,7 +946,6 @@ void computeProductPotentials(Vector3* force,
 	    SimplePotential& p = potentialList[ productPotential_list[i].x + j ];
 	    if (tmp_force != 0) {
 		// TODO add energy
-		// TODO make it work with replicas
 		p.apply_force(pos,sys, force, &productPotentialParticles[part_idx], tmp_force);
 	    }
 	    part_idx += p.type==BOND? 2: p.type==ANGLE? 3: 4;
-- 
GitLab