diff --git a/src/Angle.cu b/src/Angle.cu
index 28facb728927e782b5cb50697fdb2a80618557aa..220118759eb45508b49bea6f02b24eb9ac820825 100644
--- a/src/Angle.cu
+++ b/src/Angle.cu
@@ -15,15 +15,15 @@ String Angle::toString()
return String("ANGLE ") + ind1 + " " + ind2 + " " + ind3 + " " + fileName;
}
-void BondAngle::print()
-{
-// printf("about to print fileName %p\n", fileName.val());
-// fileName.print();
- printf("BONDANGLE (%d %d %d) %s; %s; %s\n", ind1, ind2, ind3, angleFileName.val(), bondFileName1.val(), bondFileName2.val());
-}
+// void BondAngle::print()
+// {
+// // printf("about to print fileName %p\n", fileName.val());
+// // fileName.print();
+// printf("BONDANGLE (%d %d %d) %s; %s; %s\n", ind1, ind2, ind3, angleFileName.val(), bondFileName1.val(), bondFileName2.val());
+// }
-String BondAngle::toString()
-{
- return String("BONDANGLE ") + ind1 + " " + ind2 + " " + ind3 + " " + angleFileName + " " + bondFileName1 + " " + bondFileName2;
-}
+// String BondAngle::toString()
+// {
+// return String("BONDANGLE ") + ind1 + " " + ind2 + " " + ind3 + " " + angleFileName + " " + bondFileName1 + " " + bondFileName2;
+// }
diff --git a/src/Angle.h b/src/Angle.h
index a8794e89c099d9b7d18d294c7065246586685ebf..e03a77914ec9a22eff5a67df20f4d79416d32a88 100644
--- a/src/Angle.h
+++ b/src/Angle.h
@@ -65,14 +65,14 @@ class BondAngle
{
public:
BondAngle() {}
- BondAngle(int ind1, int ind2, int ind3, String angleFileName, String bondFileName1, String bondFileName2) :
- ind1(ind1), ind2(ind2), ind3(ind3), angleFileName(angleFileName), bondFileName1(bondFileName1), bondFileName2(bondFileName2), tabFileIndex1(-1), tabFileIndex2(-1), tabFileIndex3(-1) { }
+ BondAngle(int ind1, int ind2, int ind3, int ind4, String angleFileName1, String bondFileName, String angleFileName2) :
+ ind1(ind1), ind2(ind2), ind3(ind3), ind4(ind4), angleFileName1(angleFileName1), bondFileName(bondFileName), angleFileName2(angleFileName2), tabFileIndex1(-1), tabFileIndex2(-1), tabFileIndex3(-1) { }
- int ind1, ind2, ind3;
+ int ind1, ind2, ind3, ind4;
- String angleFileName;
- String bondFileName1;
- String bondFileName2;
+ String angleFileName1;
+ String bondFileName;
+ String angleFileName2;
// tabFileIndex will be assigned after ComputeForce loads the
// TabulatedAnglePotentials. The tabefileIndex is used by ComputeForce to
// discern which TabulatedAnglePotential this Angle uses.
@@ -80,38 +80,35 @@ public:
int tabFileIndex2;
int tabFileIndex3;
- inline BondAngle(const BondAngle& a) : ind1(a.ind1), ind2(a.ind2), ind3(a.ind3),
- angleFileName(a.angleFileName),
- bondFileName1(a.bondFileName1),
- bondFileName2(a.bondFileName2),
- tabFileIndex1(a.tabFileIndex1),
- tabFileIndex2(a.tabFileIndex2),
- tabFileIndex3(a.tabFileIndex3) { }
-
- HOST DEVICE inline float calcAngle(Vector3* pos, BaseGrid* sys) {
- const Vector3& posa = pos[ind1];
- const Vector3& posb = pos[ind2];
- const Vector3& posc = pos[ind3];
- const float distab = sys->wrapDiff(posa - posb).length();
- const float distbc = sys->wrapDiff(posb - posc).length();
- const float distac = sys->wrapDiff(posc - posa).length();
- float cos = (distbc * distbc + distab * distab - distac * distac)
- / (2.0f * distbc * distab);
- if (cos < -1.0f) cos = -1.0f;
- else if (cos > 1.0f) cos = 1.0f;
- float angle = acos(cos);
- return angle;
- }
-
- HOST DEVICE inline int getIndex(int index) {
- if (index == ind1) return 1;
- if (index == ind2) return 2;
- if (index == ind3) return 3;
- return -1;
- }
-
- String toString();
- void print();
+ inline BondAngle(const BondAngle& a) : ind1(a.ind1), ind2(a.ind2), ind3(a.ind3), ind4(a.ind4),
+ angleFileName1(a.angleFileName1), bondFileName(a.bondFileName), angleFileName2(a.angleFileName2),
+ tabFileIndex1(a.tabFileIndex1), tabFileIndex2(a.tabFileIndex2), tabFileIndex3(a.tabFileIndex3) { }
+
+ // HOST DEVICE inline float calcAngle(Vector3* pos, BaseGrid* sys) {
+ // const Vector3& posa = pos[ind1];
+ // const Vector3& posb = pos[ind2];
+ // const Vector3& posc = pos[ind3];
+ // const float distab = sys->wrapDiff(posa - posb).length();
+ // const float distbc = sys->wrapDiff(posb - posc).length();
+ // const float distac = sys->wrapDiff(posc - posa).length();
+ // float cos = (distbc * distbc + distab * distab - distac * distac)
+ // / (2.0f * distbc * distab);
+ // if (cos < -1.0f) cos = -1.0f;
+ // else if (cos > 1.0f) cos = 1.0f;
+ // float angle = acos(cos);
+ // return angle;
+ // }
+
+ // HOST DEVICE inline int getIndex(int index) {
+ // if (index == ind1) return 1;
+ // if (index == ind2) return 2;
+ // if (index == ind3) return 3;
+ // if (index == ind4) return 4;
+ // return -1;
+ // }
+
+ // String toString();
+ // void print();
};
#endif
diff --git a/src/ComputeForce.cu b/src/ComputeForce.cu
index c66fb14b3bd4f76ddc2c5651485f3bd4d0064bf0..60bc945408cca1f2eaf71cbd9c7d6dec68460c07 100644
--- a/src/ComputeForce.cu
+++ b/src/ComputeForce.cu
@@ -419,10 +419,8 @@ bool ComputeForce::addBondPotential(String fileName, int ind, Bond bonds[], Bond
for (int i = 0; i < numBondAngles; i++)
{
- if (bondAngles[i].bondFileName1 == fileName)
+ if (bondAngles[i].bondFileName == fileName)
bondAngles[i].tabFileIndex2 = ind;
- if (bondAngles[i].bondFileName2 == fileName)
- bondAngles[i].tabFileIndex3 = ind;
}
gpuErrchk(cudaMemcpyAsync(bonds_d, bonds, sizeof(Bond) * numBonds, cudaMemcpyHostToDevice));
@@ -462,10 +460,12 @@ bool ComputeForce::addAnglePotential(String fileName, int ind, Angle* angles, Bo
if (angles[i].fileName == fileName)
angles[i].tabFileIndex = ind;
- for (int i = 0; i < numBondAngles; i++)
- if (bondAngles[i].angleFileName == fileName)
- bondAngles[i].tabFileIndex1 = ind;
-
+ for (int i = 0; i < numBondAngles; i++) {
+ if (bondAngles[i].angleFileName1 == fileName)
+ bondAngles[i].tabFileIndex1 = ind;
+ if (bondAngles[i].angleFileName2 == fileName)
+ bondAngles[i].tabFileIndex3 = ind;
+ }
gpuErrchk(cudaMemcpy(angles_d, angles, sizeof(Angle) * numAngles,
cudaMemcpyHostToDevice));
return true;
@@ -979,7 +979,7 @@ void ComputeForce::copyToCUDA(int simNum, int *type, Bond* bonds, int2* bondMap,
// }
// }
-void ComputeForce::copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *dihedralList, int *dihedralPotList, int2* bondAngleList) {
+void ComputeForce::copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *dihedralList, int *dihedralPotList, int4* bondAngleList) {
size_t size;
@@ -1007,7 +1007,7 @@ void ComputeForce::copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *d
}
if (numBondAngles > 0) {
- size = 3*numBondAngles * numReplicas * sizeof(int2);
+ size = 2*numBondAngles * numReplicas * sizeof(int4);
gpuErrchk( cudaMalloc( &bondAngleList_d, size ) );
gpuErrchk( cudaMemcpyAsync( bondAngleList_d, bondAngleList, size, cudaMemcpyHostToDevice) );
}
diff --git a/src/ComputeForce.cuh b/src/ComputeForce.cuh
index 16de7427b2678600fe554002217cfe4503598317..bbbce8c38d5617e96ec3ceaf63eba7b0659cbfa7 100644
--- a/src/ComputeForce.cuh
+++ b/src/ComputeForce.cuh
@@ -869,20 +869,21 @@ __global__
void computeTabulatedBondAngles(Vector3* force,
Vector3* __restrict__ pos,
BaseGrid* __restrict__ sys,
- int numBondAngles, int2* __restrict__ bondAngleList_d, TabulatedAnglePotential** tableAngle,
+ int numBondAngles, int4* __restrict__ bondAngleList_d, TabulatedAnglePotential** tableAngle,
TabulatedPotential** tableBond,
float* energy, bool get_energy) {
// Loop over ALL angles in ALL replicas
for (int i = threadIdx.x+blockIdx.x*blockDim.x; i<numBondAngles; i+=blockDim.x*gridDim.x) {
- int atom1 = bondAngleList_d[3*i ].x;
- int atom2 = bondAngleList_d[3*i ].y;
- int atom3 = bondAngleList_d[3*i+1].x;
+ int atom1 = bondAngleList_d[2*i].x;
+ int atom2 = bondAngleList_d[2*i].y;
+ int atom3 = bondAngleList_d[2*i].z;
+ int atom4 = bondAngleList_d[2*i].w;
- int angleInd = bondAngleList_d[3*i+1].y;
- int bondInd1 = bondAngleList_d[3*i+2].x;
- int bondInd2 = bondAngleList_d[3*i+2].y;
+ int angleInd1 = bondAngleList_d[2*i+1].x;
+ int bondInd = bondAngleList_d[2*i+1].y;
+ int angleInd2 = bondAngleList_d[2*i+1].z;
- computeBondAngle(tableAngle[ angleInd ], tableBond[ bondInd1 ], tableBond[ bondInd2 ], sys, force, pos, atom1, atom2, atom3, energy, get_energy);
+ computeBondAngle(tableAngle[ angleInd1 ], tableBond[ bondInd ], tableAngle[ angleInd2 ], sys, force, pos, atom1, atom2, atom3, atom4, energy, get_energy);
}
}
diff --git a/src/ComputeForce.h b/src/ComputeForce.h
index 08269d511d7aec24af27700c092dc484fbab7845..d4b83228b796b4b5d4d3fcead66e8dfe0d375039 100644
--- a/src/ComputeForce.h
+++ b/src/ComputeForce.h
@@ -80,7 +80,7 @@ public:
void copyToCUDA(Vector3* forceInternal, Vector3* pos, Vector3* mom, float* random);
// void createBondList(int3 *bondList);
- void copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *dihedralList, int *dihedralPotList, int2 *bondAngleList);
+ void copyBondedListsToGPU(int3 *bondList, int4 *angleList, int4 *dihedralList, int *dihedralPotList, int4 *bondAngleList);
//MLog: because of the move of a lot of private variables, some functions get starved necessary memory access to these variables, below is a list of functions that return the specified private variable.
std::vector<Vector3*> getPos_d()
@@ -251,7 +251,7 @@ private:
int numBondAngles;
BondAngle* bondAngles_d;
- int2* bondAngleList_d;
+ int4* bondAngleList_d;
int3* bondList_d;
int4* angleList_d;
diff --git a/src/Configuration.cpp b/src/Configuration.cpp
index a6f06b8e9cc1078a9b3492287420903391d94783..e12f81153fb9f753b058100ab860f8bd98f0bd41 100644
--- a/src/Configuration.cpp
+++ b/src/Configuration.cpp
@@ -1962,7 +1962,7 @@ void Configuration::readBondAngles() {
// Legitimate BONDANGLE inputs have 7 tokens
// BONDANGLE | INDEX1 | INDEX2 | INDEX3 | ANGLE_FILENAME | BOND_FILENAME1 | BOND_FILENAME2
// Any angle input line without exactly 5 tokens should be discarded
- if (numTokens != 7) {
+ if (numTokens != 8) {
printf("WARNING: Invalid bond_angle input line: %s\n", line);
continue;
}
@@ -1974,11 +1974,12 @@ void Configuration::readBondAngles() {
int ind1 = atoi(tokenList[1].val());
int ind2 = atoi(tokenList[2].val());
int ind3 = atoi(tokenList[3].val());
- String file_name1 = tokenList[4];
- String file_name2 = tokenList[5];
- String file_name3 = tokenList[6];
+ int ind4 = atoi(tokenList[4].val());
+ String file_name1 = tokenList[5];
+ String file_name2 = tokenList[6];
+ String file_name3 = tokenList[7];
//printf("file_name %s\n", file_name.val());
- if (ind1 >= num or ind2 >= num or ind3 >= num)
+ if (ind1 >= num or ind2 >= num or ind3 >= num or ind4 >= num)
continue;
if (numBondAngles >= capacity) {
@@ -1990,7 +1991,7 @@ void Configuration::readBondAngles() {
delete[] temp;
}
- BondAngle a(ind1, ind2, ind3, file_name1, file_name2, file_name3);
+ BondAngle a(ind1, ind2, ind3, ind4, file_name1, file_name2, file_name3);
bondAngles[numBondAngles++] = a;
delete[] tokenList;
}
@@ -2607,5 +2608,8 @@ bool Configuration::compare::operator()(const BondAngle& lhs, const BondAngle& r
diff = lhs.ind2 - rhs.ind2;
if (diff != 0)
return lhs.ind2 < rhs.ind2;
- return lhs.ind3 < rhs.ind3;
+ diff = lhs.ind3 - rhs.ind3;
+ if (diff != 0)
+ return lhs.ind3 < rhs.ind3;
+ return lhs.ind4 < rhs.ind4;
}
diff --git a/src/GrandBrownTown.cu b/src/GrandBrownTown.cu
index 7ea4cd27f498458cbd1e2be93bd76cd0a961ad6b..6353b8a92095b30a0f13438eb2db70283907d8b0 100644
--- a/src/GrandBrownTown.cu
+++ b/src/GrandBrownTown.cu
@@ -379,25 +379,25 @@ GrandBrownTown::GrandBrownTown(const Configuration& c, const char* outArg,
}
if (numBondAngles > 0) {
- bondAngleList = new int2[ (numBondAngles*3) * numReplicas ];
+ bondAngleList = new int4[ (numBondAngles*2) * numReplicas ];
for(int k = 0 ; k < numReplicas; k++) {
for(int i = 0; i < numBondAngles; ++i) {
if (bondAngles[i].tabFileIndex1 == -1) {
- fprintf(stderr,"Error: bondanglefile '%s' was not read with tabulatedAngleFile command.\n", bondAngles[i].angleFileName.val());
+ fprintf(stderr,"Error: bondanglefile '%s' was not read with tabulatedAngleFile command.\n", bondAngles[i].angleFileName1.val());
exit(1);
}
if (bondAngles[i].tabFileIndex2 == -1) {
- fprintf(stderr,"Error: bondanglefile1 '%s' was not read with tabulatedBondFile command.\n", bondAngles[i].bondFileName1.val());
+ fprintf(stderr,"Error: bondanglefile1 '%s' was not read with tabulatedBondFile command.\n", bondAngles[i].bondFileName.val());
exit(1);
}
if (bondAngles[i].tabFileIndex3 == -1) {
- fprintf(stderr,"Error: bondanglefile2 '%s' was not read with tabulatedBondFile command.\n", bondAngles[i].bondFileName2.val());
+ fprintf(stderr,"Error: bondanglefile2 '%s' was not read with tabulatedBondFile command.\n", bondAngles[i].angleFileName2.val());
exit(1);
}
int idx = i+k*numBondAngles;
- bondAngleList[idx*3] = make_int2( bondAngles[i].ind1+k*num, bondAngles[i].ind2+k*num );
- bondAngleList[idx*3+1] = make_int2( bondAngles[i].ind3+k*num, bondAngles[i].tabFileIndex1 );
- bondAngleList[idx*3+2] = make_int2( bondAngles[i].tabFileIndex2, bondAngles[i].tabFileIndex3 );
+ bondAngleList[idx*2] = make_int4( bondAngles[i].ind1+k*num, bondAngles[i].ind2+k*num,
+ bondAngles[i].ind3+k*num, bondAngles[i].ind4+k*num );
+ bondAngleList[idx*2+1] = make_int4( bondAngles[i].tabFileIndex1, bondAngles[i].tabFileIndex2, bondAngles[i].tabFileIndex3, -1 );
}
}
}
diff --git a/src/GrandBrownTown.h b/src/GrandBrownTown.h
index c74ca2b59f243a8ff3867e3581de46f338af2b64..ce7f19c29534db17e8b2ef8656f866b701455c67 100644
--- a/src/GrandBrownTown.h
+++ b/src/GrandBrownTown.h
@@ -277,7 +277,7 @@ private:
int numBondAngles;
BondAngle* bondAngles;
- int2 *bondAngleList;
+ int4 *bondAngleList;
//Han-Yi Chou
String particle_dynamic;
diff --git a/src/TabulatedMethods.cuh b/src/TabulatedMethods.cuh
index 7965c7736c793bbd7b0cfcf7d6229d8734247f4b..c775a6b00acd336a95615a496244b45eee9e3a36 100644
--- a/src/TabulatedMethods.cuh
+++ b/src/TabulatedMethods.cuh
@@ -2,6 +2,14 @@
// Defined elsewhere: constexpr float BD_PI = 3.1415927f;
+struct AngleForce {
+ __host__ __device__
+ AngleForce(Vector3 f1, Vector3 f3, float e) : f1(f1), f3(f3), e(e) { }
+ Vector3 f1;
+ Vector3 f3;
+ float e;
+};
+
__device__ inline void computeAngle(const TabulatedAnglePotential* __restrict__ a, const BaseGrid* __restrict__ sys, Vector3* force, const Vector3* __restrict__ pos,
const int& i, const int& j, const int& k, float* energy, bool get_energy) {
@@ -77,43 +85,31 @@ __device__ inline void computeAngle(const TabulatedAnglePotential* __restrict__
atomicAdd( &force[k], force3 );
}
-__device__ inline void computeBondAngle(const TabulatedAnglePotential* __restrict__ a,
- const TabulatedPotential* __restrict__ b1, const TabulatedPotential* __restrict__ b2,
- const BaseGrid* __restrict__ sys, Vector3* force, const Vector3* __restrict__ pos,
- const int& i, const int& j, const int& k, float* energy, bool get_energy) {
-
- // Particle's type and position
- Vector3 posa = pos[i];
- Vector3 posb = pos[j];
- Vector3 posc = pos[k];
-
- // The vectors between each pair of particles
- const Vector3 ab = sys->wrapDiff(posa - posb);
- const Vector3 bc = sys->wrapDiff(posb - posc);
- const Vector3 ac = sys->wrapDiff(posc - posa);
-
+__device__ inline AngleForce calcAngle(const TabulatedAnglePotential* __restrict__ a, const Vector3 ab, const Vector3 bc, const Vector3 ac) {
+ // // The vectors between each pair of particles
+ // const Vector3 ab = sys->wrapDiff(posa - posb);
+ // const Vector3 bc = sys->wrapDiff(posb - posc);
+ // const Vector3 ac = sys->wrapDiff(posc - posa);
+
// Find the distance between each pair of particles
float distab = ab.length2();
float distbc = bc.length2();
- EnergyForce fe_b1 = b1->compute(ab,distab);
- EnergyForce fe_b2 = b2->compute(bc,distbc);
-
const float distac2 = ac.length2();
-
- // Find the cosine of the angle we want - <ABC
+
+ // Find the cosine of the angle we want - <ABC
float cos = (distab + distbc - distac2);
distab = 1.0f/sqrt(distab); //TODO: test other functiosn
distbc = 1.0f/sqrt(distbc);
cos *= 0.5f * distbc * distab;
-
+
// If the cosine is illegitimate, set it to 1 or -1 so that acos won't fail
if (cos < -1.0f) cos = -1.0f;
if (cos > 1.0f) cos = 1.0f;
// Find the sine while we're at it.
- // Now we can use the cosine to find the actual angle (in radians)
+ // Now we can use the cosine to find the actual angle (in radians)
float angle = acos(cos);
// transform angle to units of tabulated array index
@@ -121,42 +117,74 @@ __device__ inline void computeBondAngle(const TabulatedAnglePotential* __restric
// tableAngle[0] stores the potential at angle_step
// tableAngle[1] stores the potential at angle_step * 2, etc.
- // 'home' is the index after which 'convertedAngle' would appear if it were stored in the table
+ // 'home' is the index after which 'convertedAngle' would appear if it were stored in the table
int home = int(floorf(angle));
- home = (home >= a->size) ? (a->size)-1 : home;
+ home = (home >= a->size) ? (a->size)-1 : home;
//assert(home >= 0);
//assert(home+1 < a->size);
// // Make angle the distance from [0,1) from the first index in the potential array index
// angle -= home;
-
- // Linearly interpolate the potential
+
+ // Linearly interpolate the potential
float U0 = a->pot[home];
float dUdx = (a->pot[(((home+1)==(a->size)) ? (a->size)-1 : home+1)] - U0) * a->angle_step_inv;
- float a_e = (dUdx * (angle-home)) + U0;
- if(get_energy)
- {
- float e = a_e * fe_b1.e * fe_b2.e * 0.3333333333f;
- atomicAdd( &energy[i], e);
- atomicAdd( &energy[j], e);
- atomicAdd( &energy[k], e);
- }
+ float e = ((dUdx * (angle-home)) + U0);
+
float sin = sqrtf(1.0f - cos*cos);
- dUdx /= abs(sin) > 1e-3 ? sin : 1e-3; // avoid singularity
+ dUdx /= abs(sin) > 1e-3 ? sin : 1e-3; // avoid singularity
// Calculate the forces
Vector3 force1 = -(dUdx*distab) * (ab * (cos*distab) + bc * distbc); // force on particle 1
Vector3 force3 = (dUdx*distbc) * (bc * (cos*distbc) + ab * distab); // force on particle 3
- force1 = (force1 * fe_b1.e - a_e * fe_b1.f) * fe_b2.e;
- force3 = (force3 * fe_b2.e + a_e * fe_b2.f) * fe_b1.e;
+ return AngleForce(force1,force3,e);
+}
- // assert( force1.length() < 10000.0f );
- // assert( force3.length() < 10000.0f );
+__device__ inline void computeBondAngle(const TabulatedAnglePotential* __restrict__ a1,
+ const TabulatedPotential* __restrict__ b, const TabulatedAnglePotential* __restrict__ a2,
+ const BaseGrid* __restrict__ sys, Vector3* force, const Vector3* __restrict__ pos,
+ const int& i, const int& j, const int& k, const int& l, float* energy, bool get_energy) {
- atomicAdd( &force[i], force1 );
- atomicAdd( &force[j], -(force1 + force3) );
- atomicAdd( &force[k], force3 );
+ // Particle's type and position
+ Vector3 posa = pos[i];
+ Vector3 posb = pos[j];
+ Vector3 posc = pos[k];
+ Vector3 posd = pos[l];
+
+ // The vectors between each pair of particles
+ const Vector3 ab = sys->wrapDiff(posb - posa);
+ const Vector3 bc = sys->wrapDiff(posc - posb);
+ const Vector3 ca = sys->wrapDiff(posc - posa);
+ AngleForce fe_a1 = calcAngle(a1, -ab,-bc,ca);
+
+ float distbc = bc.length2();
+ EnergyForce fe_b = b->compute(bc,distbc);
+
+ const Vector3 cd = sys->wrapDiff(posd - posc);
+ const Vector3 db = sys->wrapDiff(posd - posb);
+ AngleForce fe_a2 = calcAngle(a2, -bc,-cd,db);
+
+ if(get_energy)
+ {
+ float e = fe_a1.e * fe_b.e * fe_a2.e * 0.25f;
+ atomicAdd( &energy[i], e);
+ atomicAdd( &energy[j], e);
+ atomicAdd( &energy[k], e);
+ atomicAdd( &energy[l], e);
+ }
+ atomicAdd( &force[i], fe_a1.f1 * fe_b.e * fe_a2.e );
+ atomicAdd( &force[j],
+ -(fe_a1.f1 + fe_a1.f3) * fe_b.e * fe_a2.e
+ + fe_b.f * fe_a1.e * fe_a2.e
+ + fe_a2.f1 * fe_b.e * fe_a1.e
+ );
+ atomicAdd( &force[k],
+ fe_a1.f3 * fe_b.e * fe_a2.e
+ - fe_b.f * fe_a1.e * fe_a2.e
+ - (fe_a2.f1 + fe_a2.f3) * fe_b.e * fe_a1.e
+ );
+ atomicAdd( &force[l], fe_a2.f3 * fe_b.e * fe_a1.e );
}