Commit 8f532e31 authored by peastman's avatar peastman
Browse files

Optimization

parent f854d108
...@@ -356,8 +356,8 @@ RealOpenMM CpuGayBerneForce::computeOneInteraction(int particle1, int particle2, ...@@ -356,8 +356,8 @@ RealOpenMM CpuGayBerneForce::computeOneInteraction(int particle1, int particle2,
// Interactions between two point particles can be computed more easily. // Interactions between two point particles can be computed more easily.
if (particles[particle1].isPointParticle && particles[particle2].isPointParticle) { if (particles[particle1].isPointParticle && particles[particle2].isPointParticle) {
RealOpenMM sig2 = sigma*rInv; RealOpenMM sig = sigma*rInv;
sig2 *= sig2; RealOpenMM sig2 = sig*sig;
RealOpenMM sig6 = sig2*sig2*sig2; RealOpenMM sig6 = sig2*sig2*sig2;
RealOpenMM energy = 4*epsilon*(sig6-1)*sig6; RealOpenMM energy = 4*epsilon*(sig6-1)*sig6;
RealVec force = drUnit*(switchValue*4*epsilon*(12*sig6 - 6)*sig6*rInv - energy*switchDeriv); RealVec force = drUnit*(switchValue*4*epsilon*(12*sig6 - 6)*sig6*rInv - energy*switchDeriv);
...@@ -422,51 +422,33 @@ RealOpenMM CpuGayBerneForce::computeOneInteraction(int particle1, int particle2, ...@@ -422,51 +422,33 @@ RealOpenMM CpuGayBerneForce::computeOneInteraction(int particle1, int particle2,
RealVec scale = RealVec(p.rx*p.rx, p.ry*p.ry, p.rz*p.rz)*(-0.5*eta/detG12); RealVec scale = RealVec(p.rx*p.rx, p.ry*p.ry, p.rz*p.rz)*(-0.5*eta/detG12);
Matrix D; Matrix D;
RealOpenMM (&d)[3][3] = D.v; RealOpenMM (&d)[3][3] = D.v;
d[0][0] = scale[0]*(g12[1][2]*g12[0][1]*a[0][2] + 2*g12[1][1]*g12[2][2]*a[0][0] - d[0][0] = scale[0]*(2*a[0][0]*(g12[1][1]*g12[2][2] - g12[1][2]*g12[2][1]) +
g12[1][1]*a[0][2]*g12[0][2] - 2*g12[1][2]*a[0][0]*g12[2][1] + a[0][2]*(g12[1][2]*g12[0][1] + g12[1][0]*g12[2][1] - g12[1][1]*(g12[0][2] + g12[2][0])) +
a[0][1]*g12[0][2]*g12[2][1] - a[0][1]*g12[0][1]*g12[2][2] - a[0][1]*(g12[0][2]*g12[2][1] + g12[2][0]*g12[1][2] - g12[2][2]*(g12[0][1] + g12[1][0])));
g12[1][0]*g12[2][2]*a[0][1] + g12[2][0]*g12[1][2]*a[0][1] + d[0][1] = scale[0]*( a[0][0]*(g12[0][2]*g12[2][1] + g12[2][0]*g12[1][2] - g12[2][2]*(g12[0][1] + g12[1][0])) +
g12[1][0]*a[0][2]*g12[2][1] - a[0][2]*g12[2][0]*g12[1][1]); 2*a[0][1]*(g12[0][0]*g12[2][2] - g12[2][0]*g12[0][2]) +
d[0][1] = scale[0]*( g12[0][2]*a[0][0]*g12[2][1] - g12[2][2]*a[0][0]*g12[0][1] + a[0][2]*(g12[1][0]*g12[0][2] + g12[2][0]*g12[0][1] - g12[0][0]*(g12[1][2] + g12[2][1])));
2*g12[0][0]*g12[2][2]*a[0][1] - g12[0][0]*a[0][2]*g12[1][2] - d[0][2] = scale[0]*( a[0][0]*(g12[0][1]*g12[1][2] + g12[1][0]*g12[2][1] - g12[1][1]*(g12[0][2] + g12[2][0])) +
2*g12[2][0]*g12[0][2]*a[0][1] + a[0][2]*g12[1][0]*g12[0][2] - a[0][1]*(g12[1][0]*g12[0][2] + g12[2][0]*g12[0][1] - g12[0][0]*(g12[1][2] + g12[2][1])) +
g12[2][2]*g12[1][0]*a[0][0] + g12[2][0]*a[0][0]*g12[1][2] + 2*a[0][2]*(g12[1][1]*g12[0][0] - g12[1][0]*g12[0][1]));
g12[2][0]*a[0][2]*g12[0][1] - a[0][2]*g12[0][0]*g12[2][1]); d[1][0] = scale[1]*(2*a[1][0]*(g12[1][1]*g12[2][2] - g12[1][2]*g12[2][1]) +
d[0][2] = scale[0]*( g12[0][1]*g12[1][2]*a[0][0] - g12[0][2]*a[0][0]*g12[1][1] - a[1][1]*(g12[0][2]*g12[2][1] + g12[2][0]*g12[1][2] - g12[2][2]*(g12[0][1] + g12[1][0])) +
g12[0][0]*g12[1][2]*a[0][1] + g12[1][0]*g12[0][2]*a[0][1] - a[1][2]*(g12[1][2]*g12[0][1] + g12[1][0]*g12[2][1] - g12[1][1]*(g12[0][2] + g12[2][0])));
a[0][1]*g12[0][0]*g12[2][1] - g12[2][0]*g12[1][1]*a[0][0] + d[1][1] = scale[1]*( a[1][0]*(g12[0][2]*g12[2][1] + g12[2][0]*g12[1][2] - g12[2][2]*(g12[0][1] + g12[1][0])) +
2*g12[1][1]*g12[0][0]*a[0][2] - 2*g12[1][0]*a[0][2]*g12[0][1] + 2*a[1][1]*(g12[2][2]*g12[0][0] - g12[2][0]*g12[0][2]) +
g12[1][0]*g12[2][1]*a[0][0] + g12[2][0]*a[0][1]*g12[0][1]); a[1][2]*(g12[1][0]*g12[0][2] + g12[0][1]*g12[2][0] - g12[0][0]*(g12[1][2] + g12[2][1])));
d[1][0] = scale[1]*(-g12[1][1]*a[1][2]*g12[0][2] + 2*g12[1][1]*g12[2][2]*a[1][0] + d[1][2] = scale[1]*( a[1][0]*(g12[0][1]*g12[1][2] + g12[1][0]*g12[2][1] - g12[1][1]*(g12[0][2] + g12[2][0])) +
g12[1][2]*g12[0][1]*a[1][2] - 2*g12[1][2]*a[1][0]*g12[2][1] + a[1][1]*(g12[1][0]*g12[0][2] + g12[0][1]*g12[2][0] - g12[0][0]*(g12[1][2] + g12[2][1])) +
a[1][1]*g12[0][2]*g12[2][1] - a[1][1]*g12[0][1]*g12[2][2] - 2*a[1][2]*(g12[1][1]*g12[0][0] - g12[1][0]*g12[0][1]));
g12[1][0]*g12[2][2]*a[1][1] + g12[2][0]*g12[1][2]*a[1][1] - d[2][0] = scale[2]*(2*a[2][0]*(g12[1][1]*g12[2][2] - g12[2][1]*g12[1][2]) +
a[1][2]*g12[2][0]*g12[1][1] + g12[1][0]*a[1][2]*g12[2][1]); a[2][1]*(g12[0][2]*g12[2][1] + g12[1][2]*g12[2][0] - g12[2][2]*(g12[0][1] + g12[1][0])) +
d[1][1] = scale[1]*( g12[0][2]*a[1][0]*g12[2][1] - g12[0][1]*g12[2][2]*a[1][0] + a[2][2]*(g12[0][1]*g12[1][2] + g12[2][1]*g12[1][0] - g12[1][1]*(g12[0][2] + g12[2][0])));
2*g12[2][2]*g12[0][0]*a[1][1] - a[1][2]*g12[0][0]*g12[1][2] - d[2][1] = scale[2]*( a[2][0]*(g12[0][2]*g12[2][1] + g12[1][2]*g12[2][0] - g12[2][2]*(g12[0][1] + g12[1][0])) +
2*g12[2][0]*a[1][1]*g12[0][2] - g12[1][0]*g12[2][2]*a[1][0] + 2*a[2][1]*(g12[0][0]*g12[2][2] - g12[0][2]*g12[2][0]) +
g12[2][0]*g12[1][2]*a[1][0] + g12[1][0]*a[1][2]*g12[0][2] - a[2][2]*(g12[1][0]*g12[0][2] + g12[0][1]*g12[2][0] - g12[0][0]*(g12[1][2] + g12[2][1])));
g12[0][0]*a[1][2]*g12[2][1] + a[1][2]*g12[0][1]*g12[2][0]); d[2][2] = scale[2]*( a[2][0]*(g12[0][1]*g12[1][2] + g12[2][1]*g12[1][0] - g12[1][1]*(g12[0][2] + g12[2][0])) +
d[1][2] = scale[1]*( g12[0][1]*g12[1][2]*a[1][0] - g12[0][2]*a[1][0]*g12[1][1] - a[2][1]*(g12[1][0]*g12[0][2] + g12[2][0]*g12[0][1] - g12[0][0]*(g12[1][2] + g12[2][1])) +
g12[0][0]*g12[1][2]*a[1][1] + g12[1][0]*g12[0][2]*a[1][1] + 2*a[2][2]*(g12[1][1]*g12[0][0] - g12[1][0]*g12[0][1]));
2*g12[1][1]*g12[0][0]*a[1][2] - g12[0][0]*a[1][1]*g12[2][1] +
g12[0][1]*g12[2][0]*a[1][1] - a[1][0]*g12[2][0]*g12[1][1] -
2*g12[1][0]*g12[0][1]*a[1][2] + g12[1][0]*a[1][0]*g12[2][1]);
d[2][0] = scale[2]*( -g12[1][1]*g12[0][2]*a[2][2] + g12[0][1]*g12[1][2]*a[2][2] +
2*g12[1][1]*a[2][0]*g12[2][2] - g12[0][1]*a[2][1]*g12[2][2] +
g12[0][2]*g12[2][1]*a[2][1] - 2*a[2][0]*g12[2][1]*g12[1][2] -
g12[1][0]*a[2][1]*g12[2][2] + g12[1][2]*g12[2][0]*a[2][1] -
g12[1][1]*g12[2][0]*a[2][2] + g12[2][1]*g12[1][0]*a[2][2]);
d[2][1] = scale[2]*( -g12[0][1]*g12[2][2]*a[2][0] + g12[0][2]*a[2][0]*g12[2][1] +
2*a[2][1]*g12[0][0]*g12[2][2] - g12[1][2]*a[2][2]*g12[0][0] -
2*a[2][1]*g12[0][2]*g12[2][0] - g12[1][0]*a[2][0]*g12[2][2] +
g12[1][0]*g12[0][2]*a[2][2] + g12[1][2]*g12[2][0]*a[2][0] -
g12[0][0]*a[2][2]*g12[2][1] + a[2][2]*g12[0][1]*g12[2][0]);
d[2][2] = scale[2]*( g12[0][1]*g12[1][2]*a[2][0] - g12[0][2]*a[2][0]*g12[1][1] -
g12[0][0]*g12[1][2]*a[2][1] + g12[1][0]*g12[0][2]*a[2][1] -
g12[1][1]*g12[2][0]*a[2][0] - g12[2][1]*a[2][1]*g12[0][0] +
2*g12[1][1]*a[2][2]*g12[0][0] + g12[2][1]*g12[1][0]*a[2][0] +
g12[2][0]*g12[0][1]*a[2][1] - 2*a[2][2]*g12[1][0]*g12[0][1]);
RealVec detadq; RealVec detadq;
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
detadq += RealVec(a[i][0], a[i][1], a[i][2]).cross(RealVec(d[i][0], d[i][1], d[i][2])); detadq += RealVec(a[i][0], a[i][1], a[i][2]).cross(RealVec(d[i][0], d[i][1], d[i][2]));
......
...@@ -276,51 +276,33 @@ RealOpenMM ReferenceGayBerneForce::computeOneInteraction(int particle1, int part ...@@ -276,51 +276,33 @@ RealOpenMM ReferenceGayBerneForce::computeOneInteraction(int particle1, int part
RealVec scale = RealVec(p.rx*p.rx, p.ry*p.ry, p.rz*p.rz)*(-0.5*eta/detG12); RealVec scale = RealVec(p.rx*p.rx, p.ry*p.ry, p.rz*p.rz)*(-0.5*eta/detG12);
Matrix D; Matrix D;
RealOpenMM (&d)[3][3] = D.v; RealOpenMM (&d)[3][3] = D.v;
d[0][0] = scale[0]*(g12[1][2]*g12[0][1]*a[0][2] + 2*g12[1][1]*g12[2][2]*a[0][0] - d[0][0] = scale[0]*(2*a[0][0]*(g12[1][1]*g12[2][2] - g12[1][2]*g12[2][1]) +
g12[1][1]*a[0][2]*g12[0][2] - 2*g12[1][2]*a[0][0]*g12[2][1] + a[0][2]*(g12[1][2]*g12[0][1] + g12[1][0]*g12[2][1] - g12[1][1]*(g12[0][2] + g12[2][0])) +
a[0][1]*g12[0][2]*g12[2][1] - a[0][1]*g12[0][1]*g12[2][2] - a[0][1]*(g12[0][2]*g12[2][1] + g12[2][0]*g12[1][2] - g12[2][2]*(g12[0][1] + g12[1][0])));
g12[1][0]*g12[2][2]*a[0][1] + g12[2][0]*g12[1][2]*a[0][1] + d[0][1] = scale[0]*( a[0][0]*(g12[0][2]*g12[2][1] + g12[2][0]*g12[1][2] - g12[2][2]*(g12[0][1] + g12[1][0])) +
g12[1][0]*a[0][2]*g12[2][1] - a[0][2]*g12[2][0]*g12[1][1]); 2*a[0][1]*(g12[0][0]*g12[2][2] - g12[2][0]*g12[0][2]) +
d[0][1] = scale[0]*( g12[0][2]*a[0][0]*g12[2][1] - g12[2][2]*a[0][0]*g12[0][1] + a[0][2]*(g12[1][0]*g12[0][2] + g12[2][0]*g12[0][1] - g12[0][0]*(g12[1][2] + g12[2][1])));
2*g12[0][0]*g12[2][2]*a[0][1] - g12[0][0]*a[0][2]*g12[1][2] - d[0][2] = scale[0]*( a[0][0]*(g12[0][1]*g12[1][2] + g12[1][0]*g12[2][1] - g12[1][1]*(g12[0][2] + g12[2][0])) +
2*g12[2][0]*g12[0][2]*a[0][1] + a[0][2]*g12[1][0]*g12[0][2] - a[0][1]*(g12[1][0]*g12[0][2] + g12[2][0]*g12[0][1] - g12[0][0]*(g12[1][2] + g12[2][1])) +
g12[2][2]*g12[1][0]*a[0][0] + g12[2][0]*a[0][0]*g12[1][2] + 2*a[0][2]*(g12[1][1]*g12[0][0] - g12[1][0]*g12[0][1]));
g12[2][0]*a[0][2]*g12[0][1] - a[0][2]*g12[0][0]*g12[2][1]); d[1][0] = scale[1]*(2*a[1][0]*(g12[1][1]*g12[2][2] - g12[1][2]*g12[2][1]) +
d[0][2] = scale[0]*( g12[0][1]*g12[1][2]*a[0][0] - g12[0][2]*a[0][0]*g12[1][1] - a[1][1]*(g12[0][2]*g12[2][1] + g12[2][0]*g12[1][2] - g12[2][2]*(g12[0][1] + g12[1][0])) +
g12[0][0]*g12[1][2]*a[0][1] + g12[1][0]*g12[0][2]*a[0][1] - a[1][2]*(g12[1][2]*g12[0][1] + g12[1][0]*g12[2][1] - g12[1][1]*(g12[0][2] + g12[2][0])));
a[0][1]*g12[0][0]*g12[2][1] - g12[2][0]*g12[1][1]*a[0][0] + d[1][1] = scale[1]*( a[1][0]*(g12[0][2]*g12[2][1] + g12[2][0]*g12[1][2] - g12[2][2]*(g12[0][1] + g12[1][0])) +
2*g12[1][1]*g12[0][0]*a[0][2] - 2*g12[1][0]*a[0][2]*g12[0][1] + 2*a[1][1]*(g12[2][2]*g12[0][0] - g12[2][0]*g12[0][2]) +
g12[1][0]*g12[2][1]*a[0][0] + g12[2][0]*a[0][1]*g12[0][1]); a[1][2]*(g12[1][0]*g12[0][2] + g12[0][1]*g12[2][0] - g12[0][0]*(g12[1][2] + g12[2][1])));
d[1][0] = scale[1]*(-g12[1][1]*a[1][2]*g12[0][2] + 2*g12[1][1]*g12[2][2]*a[1][0] + d[1][2] = scale[1]*( a[1][0]*(g12[0][1]*g12[1][2] + g12[1][0]*g12[2][1] - g12[1][1]*(g12[0][2] + g12[2][0])) +
g12[1][2]*g12[0][1]*a[1][2] - 2*g12[1][2]*a[1][0]*g12[2][1] + a[1][1]*(g12[1][0]*g12[0][2] + g12[0][1]*g12[2][0] - g12[0][0]*(g12[1][2] + g12[2][1])) +
a[1][1]*g12[0][2]*g12[2][1] - a[1][1]*g12[0][1]*g12[2][2] - 2*a[1][2]*(g12[1][1]*g12[0][0] - g12[1][0]*g12[0][1]));
g12[1][0]*g12[2][2]*a[1][1] + g12[2][0]*g12[1][2]*a[1][1] - d[2][0] = scale[2]*(2*a[2][0]*(g12[1][1]*g12[2][2] - g12[2][1]*g12[1][2]) +
a[1][2]*g12[2][0]*g12[1][1] + g12[1][0]*a[1][2]*g12[2][1]); a[2][1]*(g12[0][2]*g12[2][1] + g12[1][2]*g12[2][0] - g12[2][2]*(g12[0][1] + g12[1][0])) +
d[1][1] = scale[1]*( g12[0][2]*a[1][0]*g12[2][1] - g12[0][1]*g12[2][2]*a[1][0] + a[2][2]*(g12[0][1]*g12[1][2] + g12[2][1]*g12[1][0] - g12[1][1]*(g12[0][2] + g12[2][0])));
2*g12[2][2]*g12[0][0]*a[1][1] - a[1][2]*g12[0][0]*g12[1][2] - d[2][1] = scale[2]*( a[2][0]*(g12[0][2]*g12[2][1] + g12[1][2]*g12[2][0] - g12[2][2]*(g12[0][1] + g12[1][0])) +
2*g12[2][0]*a[1][1]*g12[0][2] - g12[1][0]*g12[2][2]*a[1][0] + 2*a[2][1]*(g12[0][0]*g12[2][2] - g12[0][2]*g12[2][0]) +
g12[2][0]*g12[1][2]*a[1][0] + g12[1][0]*a[1][2]*g12[0][2] - a[2][2]*(g12[1][0]*g12[0][2] + g12[0][1]*g12[2][0] - g12[0][0]*(g12[1][2] + g12[2][1])));
g12[0][0]*a[1][2]*g12[2][1] + a[1][2]*g12[0][1]*g12[2][0]); d[2][2] = scale[2]*( a[2][0]*(g12[0][1]*g12[1][2] + g12[2][1]*g12[1][0] - g12[1][1]*(g12[0][2] + g12[2][0])) +
d[1][2] = scale[1]*( g12[0][1]*g12[1][2]*a[1][0] - g12[0][2]*a[1][0]*g12[1][1] - a[2][1]*(g12[1][0]*g12[0][2] + g12[2][0]*g12[0][1] - g12[0][0]*(g12[1][2] + g12[2][1])) +
g12[0][0]*g12[1][2]*a[1][1] + g12[1][0]*g12[0][2]*a[1][1] + 2*a[2][2]*(g12[1][1]*g12[0][0] - g12[1][0]*g12[0][1]));
2*g12[1][1]*g12[0][0]*a[1][2] - g12[0][0]*a[1][1]*g12[2][1] +
g12[0][1]*g12[2][0]*a[1][1] - a[1][0]*g12[2][0]*g12[1][1] -
2*g12[1][0]*g12[0][1]*a[1][2] + g12[1][0]*a[1][0]*g12[2][1]);
d[2][0] = scale[2]*( -g12[1][1]*g12[0][2]*a[2][2] + g12[0][1]*g12[1][2]*a[2][2] +
2*g12[1][1]*a[2][0]*g12[2][2] - g12[0][1]*a[2][1]*g12[2][2] +
g12[0][2]*g12[2][1]*a[2][1] - 2*a[2][0]*g12[2][1]*g12[1][2] -
g12[1][0]*a[2][1]*g12[2][2] + g12[1][2]*g12[2][0]*a[2][1] -
g12[1][1]*g12[2][0]*a[2][2] + g12[2][1]*g12[1][0]*a[2][2]);
d[2][1] = scale[2]*( -g12[0][1]*g12[2][2]*a[2][0] + g12[0][2]*a[2][0]*g12[2][1] +
2*a[2][1]*g12[0][0]*g12[2][2] - g12[1][2]*a[2][2]*g12[0][0] -
2*a[2][1]*g12[0][2]*g12[2][0] - g12[1][0]*a[2][0]*g12[2][2] +
g12[1][0]*g12[0][2]*a[2][2] + g12[1][2]*g12[2][0]*a[2][0] -
g12[0][0]*a[2][2]*g12[2][1] + a[2][2]*g12[0][1]*g12[2][0]);
d[2][2] = scale[2]*( g12[0][1]*g12[1][2]*a[2][0] - g12[0][2]*a[2][0]*g12[1][1] -
g12[0][0]*g12[1][2]*a[2][1] + g12[1][0]*g12[0][2]*a[2][1] -
g12[1][1]*g12[2][0]*a[2][0] - g12[2][1]*a[2][1]*g12[0][0] +
2*g12[1][1]*a[2][2]*g12[0][0] + g12[2][1]*g12[1][0]*a[2][0] +
g12[2][0]*g12[0][1]*a[2][1] - 2*a[2][2]*g12[1][0]*g12[0][1]);
RealVec detadq; RealVec detadq;
for (int i = 0; i < 3; i++) for (int i = 0; i < 3; i++)
detadq += RealVec(a[i][0], a[i][1], a[i][2]).cross(RealVec(d[i][0], d[i][1], d[i][2])); detadq += RealVec(a[i][0], a[i][1], a[i][2]).cross(RealVec(d[i][0], d[i][1], d[i][2]));
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment