Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Jump to bottom

Optimised pair_lj_charmm_coul* code #124

Merged
merged 3 commits into from Jul 27, 2016
Merged

Optimised pair_lj_charmm_coul* code #124

Show file tree
Hide file tree
Changes from all commits
Commits
Show all changes
3 commits
Select commit Hold shift + click to select a range
b928f98
Added optimised pair_lj_* code from master
Jul 26, 2016
9670f5d
Remove the long_14* files, not relevant for this feature branch
Jul 27, 2016
46f034d
Added back in denom_lj variable, which is used by some subclasses in …
Jul 27, 2016
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Jump to
Jump to file
Failed to load files.
Diff view
Diff view
79 changes: 33 additions & 46 deletions src/KSPACE/pair_lj_charmm_coul_long.cpp
View file
Open in desktop
Expand Up @@ -172,9 +172,9 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
Expand Down Expand Up @@ -206,7 +206,7 @@ void PairLJCharmmCoulLong::compute(int eflag, int vflag)
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
evdwl *= switch1;
}
evdwl *= factor_lj;
Expand Down Expand Up @@ -247,13 +247,6 @@ void PairLJCharmmCoulLong::compute_inner()
numneigh = listinner->numneigh;
firstneigh = listinner->firstneigh;

double cut_out_on = cut_respa[0];
double cut_out_off = cut_respa[1];

double cut_out_diff = cut_out_off - cut_out_on;
double cut_out_on_sq = cut_out_on*cut_out_on;
double cut_out_off_sq = cut_out_off*cut_out_off;

// loop over neighbors of my atoms

for (ii = 0; ii < inum; ii++) {
Expand Down Expand Up @@ -289,7 +282,7 @@ void PairLJCharmmCoulLong::compute_inner()
fpair = (forcecoul + factor_lj*forcelj) * r2inv;

if (rsq > cut_out_on_sq) {
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
rsw = (sqrt(rsq) - cut_out_on)*cut_out_diff_inv;
fpair *= 1.0 + rsw*rsw*(2.0*rsw-3.0);
}

Expand Down Expand Up @@ -332,18 +325,6 @@ void PairLJCharmmCoulLong::compute_middle()
numneigh = listmiddle->numneigh;
firstneigh = listmiddle->firstneigh;

double cut_in_off = cut_respa[0];
double cut_in_on = cut_respa[1];
double cut_out_on = cut_respa[2];
double cut_out_off = cut_respa[3];

double cut_in_diff = cut_in_on - cut_in_off;
double cut_out_diff = cut_out_off - cut_out_on;
double cut_in_off_sq = cut_in_off*cut_in_off;
double cut_in_on_sq = cut_in_on*cut_in_on;
double cut_out_on_sq = cut_out_on*cut_out_on;
double cut_out_off_sq = cut_out_off*cut_out_off;

// loop over neighbors of my atoms

for (ii = 0; ii < inum; ii++) {
Expand Down Expand Up @@ -378,20 +359,20 @@ void PairLJCharmmCoulLong::compute_middle()
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}

fpair = (forcecoul + factor_lj*forcelj) * r2inv;
if (rsq < cut_in_on_sq) {
rsw = (sqrt(rsq) - cut_in_off)/cut_in_diff;
rsw = (sqrt(rsq) - cut_in_off)*cut_in_diff_inv;
fpair *= rsw*rsw*(3.0 - 2.0*rsw);
}
if (rsq > cut_out_on_sq) {
rsw = (sqrt(rsq) - cut_out_on)/cut_out_diff;
rsw = (sqrt(rsq) - cut_out_on)*cut_out_diff_inv;
fpair *= 1.0 + rsw*rsw*(2.0*rsw - 3.0);
}

Expand Down Expand Up @@ -441,13 +422,6 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
numneigh = listouter->numneigh;
firstneigh = listouter->firstneigh;

double cut_in_off = cut_respa[2];
double cut_in_on = cut_respa[3];

double cut_in_diff = cut_in_on - cut_in_off;
double cut_in_off_sq = cut_in_off*cut_in_off;
double cut_in_on_sq = cut_in_on*cut_in_on;

// loop over neighbors of my atoms

for (ii = 0; ii < inum; ii++) {
Expand Down Expand Up @@ -518,9 +492,9 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
Expand Down Expand Up @@ -562,7 +536,7 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
evdwl *= switch1;
}
evdwl *= factor_lj;
Expand Down Expand Up @@ -590,9 +564,9 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
Expand All @@ -601,9 +575,9 @@ void PairLJCharmmCoulLong::compute_outer(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
Expand Down Expand Up @@ -759,14 +733,27 @@ void PairLJCharmmCoulLong::init_style()
cut_coulsq = cut_coul * cut_coul;
cut_bothsq = MAX(cut_ljsq,cut_coulsq);

denom_lj = (cut_ljsq-cut_lj_innersq) * (cut_ljsq-cut_lj_innersq) *
(cut_ljsq-cut_lj_innersq);
denom_lj_inv = 1.0 / ( (cut_ljsq-cut_lj_innersq) * (cut_ljsq-cut_lj_innersq) *
(cut_ljsq-cut_lj_innersq) );

// set & error check interior rRESPA cutoffs

if (strstr(update->integrate_style,"respa") &&
((Respa *) update->integrate)->level_inner >= 0) {
cut_respa = ((Respa *) update->integrate)->cutoff;
cut_in_off = cut_respa[0];
cut_in_on = cut_respa[1];
cut_out_on = cut_respa[2];
cut_out_off = cut_respa[3];

cut_in_diff = cut_in_on - cut_in_off;
cut_out_diff = cut_out_off - cut_out_on;
cut_in_diff_inv = 1.0 / (cut_in_diff);
cut_out_diff_inv = 1.0 / (cut_out_diff);
cut_in_off_sq = cut_in_off*cut_in_off;
cut_in_on_sq = cut_in_on*cut_in_on;
cut_out_on_sq = cut_out_on*cut_out_on;
cut_out_off_sq = cut_out_off*cut_out_off;
if (MIN(cut_lj,cut_coul) < cut_respa[3])
error->all(FLERR,"Pair cutoff < Respa interior cutoff");
if (cut_lj_inner < cut_respa[1])
Expand Down Expand Up @@ -992,9 +979,9 @@ double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype,
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
Expand All @@ -1017,7 +1004,7 @@ double PairLJCharmmCoulLong::single(int i, int j, int itype, int jtype,
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
philj *= switch1;
}
eng += factor_lj*philj;
Expand Down
6 changes: 5 additions & 1 deletion src/KSPACE/pair_lj_charmm_coul_long.h
View file
Open in desktop
Expand Up @@ -54,7 +54,11 @@ class PairLJCharmmCoulLong : public Pair {
double cut_lj_innersq,cut_ljsq;
double cut_coul,cut_coulsq;
double cut_bothsq;
double denom_lj;
double cut_in_off, cut_in_on, cut_out_off, cut_out_on;
double cut_in_diff, cut_out_diff;
double cut_in_diff_inv, cut_out_diff_inv;
double cut_in_off_sq, cut_in_on_sq, cut_out_off_sq, cut_out_on_sq;
double denom_lj, denom_lj_inv;
double **epsilon,**sigma,**eps14,**sigma14;
double **lj1,**lj2,**lj3,**lj4,**offset;
double **lj14_1,**lj14_2,**lj14_3,**lj14_4;
Expand Down
26 changes: 13 additions & 13 deletions src/KSPACE/pair_lj_charmm_coul_msm.cpp
View file
Open in desktop
Expand Up @@ -161,9 +161,9 @@ void PairLJCharmmCoulMSM::compute(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
Expand Down Expand Up @@ -221,7 +221,7 @@ void PairLJCharmmCoulMSM::compute(int eflag, int vflag)
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
evdwl *= switch1;
}
evdwl *= factor_lj;
Expand Down Expand Up @@ -358,9 +358,9 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
Expand Down Expand Up @@ -402,7 +402,7 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag)
evdwl = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
evdwl *= switch1;
}
evdwl *= factor_lj;
Expand Down Expand Up @@ -430,9 +430,9 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
Expand All @@ -441,9 +441,9 @@ void PairLJCharmmCoulMSM::compute_outer(int eflag, int vflag)
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
Expand Down Expand Up @@ -499,9 +499,9 @@ double PairLJCharmmCoulMSM::single(int i, int j, int itype, int jtype,
forcelj = r6inv * (lj1[itype][jtype]*r6inv - lj2[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
switch2 = 12.0*rsq * (cut_ljsq-rsq) *
(rsq-cut_lj_innersq) / denom_lj;
(rsq-cut_lj_innersq) * denom_lj_inv;
philj = r6inv * (lj3[itype][jtype]*r6inv - lj4[itype][jtype]);
forcelj = forcelj*switch1 + philj*switch2;
}
Expand All @@ -524,7 +524,7 @@ double PairLJCharmmCoulMSM::single(int i, int j, int itype, int jtype,
philj = r6inv*(lj3[itype][jtype]*r6inv-lj4[itype][jtype]);
if (rsq > cut_lj_innersq) {
switch1 = (cut_ljsq-rsq) * (cut_ljsq-rsq) *
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) / denom_lj;
(cut_ljsq + 2.0*rsq - 3.0*cut_lj_innersq) * denom_lj_inv;
philj *= switch1;
}
eng += factor_lj*philj;
Expand Down
14 changes: 6 additions & 8 deletions src/USER-MISC/pair_lj_sf_dipole_sf.cpp
View file
Open in desktop
Expand Up @@ -133,8 +133,9 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag)

if (rsq < cut_coulsq[itype][jtype]) {

rcutcoul2inv=1.0/cut_coulsq[itype][jtype];
if (qtmp != 0.0 && q[j] != 0.0) {
pre1 = qtmp*q[j]*rinv*(r2inv-1.0/cut_coulsq[itype][jtype]);
pre1 = qtmp*q[j]*rinv*(r2inv-rcutcoul2inv);

forcecoulx += pre1*delx;
forcecouly += pre1*dely;
Expand All @@ -144,7 +145,6 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag)
if (mu[i][3] > 0.0 && mu[j][3] > 0.0) {
r3inv = r2inv*rinv;
r5inv = r3inv*r2inv;
rcutcoul2inv=1.0/cut_coulsq[itype][jtype];

pdotp = mu[i][0]*mu[j][0] + mu[i][1]*mu[j][1] + mu[i][2]*mu[j][2];
pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;
Expand All @@ -156,7 +156,7 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag)
aforcecouly = pre1*dely;
aforcecoulz = pre1*delz;

bfac = 1.0 - 4.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv) +
bfac = 1.0 - 4.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv +
3.0*rsq*rsq*rcutcoul2inv*rcutcoul2inv;
presf = 2.0 * r2inv * pidotr * pjdotr;
bforcecoulx = bfac * (pjdotr*mu[i][0]+pidotr*mu[j][0]-presf*delx);
Expand Down Expand Up @@ -187,10 +187,9 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag)
r3inv = r2inv*rinv;
r5inv = r3inv*r2inv;
pidotr = mu[i][0]*delx + mu[i][1]*dely + mu[i][2]*delz;
rcutcoul2inv=1.0/cut_coulsq[itype][jtype];
pre1 = 3.0 * q[j] * r5inv * pidotr * (1-rsq*rcutcoul2inv);
pqfac = 1.0 - 3.0*rsq*rcutcoul2inv +
2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv);
2.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv;
pre2 = q[j] * r3inv * pqfac;

forcecoulx += pre2*mu[i][0] - pre1*delx;
Expand All @@ -205,10 +204,9 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag)
r3inv = r2inv*rinv;
r5inv = r3inv*r2inv;
pjdotr = mu[j][0]*delx + mu[j][1]*dely + mu[j][2]*delz;
rcutcoul2inv=1.0/cut_coulsq[itype][jtype];
pre1 = 3.0 * qtmp * r5inv * pjdotr * (1-rsq*rcutcoul2inv);
qpfac = 1.0 - 3.0*rsq*rcutcoul2inv +
2.0*rsq*sqrt(rsq)*rcutcoul2inv*sqrt(rcutcoul2inv);
2.0*rsq*sqrt(rsq*rcutcoul2inv)*rcutcoul2inv;
pre2 = qtmp * r3inv * qpfac;

forcecoulx += pre1*delx - pre2*mu[j][0];
Expand Down Expand Up @@ -261,7 +259,7 @@ void PairLJSFDipoleSF::compute(int eflag, int vflag)

if (eflag) {
if (rsq < cut_coulsq[itype][jtype]) {
ecoul = (1.0-sqrt(rsq)/sqrt(cut_coulsq[itype][jtype]));
ecoul = (1.0-sqrt(rsq/cut_coulsq[itype][jtype]));
ecoul *= ecoul;
ecoul *= qtmp * q[j] * rinv;
if (mu[i][3] > 0.0 && mu[j][3] > 0.0)
Expand Down