Unverified Commit 28483c08 authored by Axel Kohlmeyer's avatar Axel Kohlmeyer Committed by GitHub
Browse files

Merge pull request #1746 from wmbrownIntel/gayberne-fix

Bug fix for gay-berne potential when mu != 1.0.
parents 2f9f4557 1962bc00
Loading
Loading
Loading
Loading
+42 B (97.6 KiB)

File changed.

No diff preview for this file type.

+1 −1
Original line number Diff line number Diff line
@@ -131,7 +131,7 @@ and

$$ \frac{ \partial \chi_{12} }{ \partial \mathbf{q}_i } = 4.0 \cdot
r^{-2} \cdot \mathbf{A}_i (- \mathbf{\iota}^T \cdot \mathbf{B}_i
\times \mathbf{\iota} ). $$
\times \mathbf{\iota} ) \cdot \mu \cdot \chi_{12}^{ ( \mu -1 ) / \mu}. $$

For the derivative of the $\eta$ term, we were unable to find a matrix
expression due to the determinant. Let $a_{mi}$ be the mth row of the
+3 −4
Original line number Diff line number Diff line
@@ -316,10 +316,9 @@ __kernel void k_gayberne(const __global numtyp4 *restrict x_,
        numtyp tempv[3];
        gpu_row_times3(iota,b1,tempv);
        gpu_cross3(tempv,iota,tchi);
        temp1 = (numtyp)-4.0*ir*ir;
        tchi[0] *= temp1;
        tchi[1] *= temp1;
        tchi[2] *= temp1;
        tchi[0] *= temp2;
        tchi[1] *= temp2;
        tchi[2] *= temp2;
      }

      numtyp temp2 = factor_lj*eta*chi;
+10 −11
Original line number Diff line number Diff line
@@ -644,10 +644,10 @@ double PairGayBerne::gayberne_analytic(const int i,const int j,double a1[3][3],
  dchi[2] = temp2*(iota[2]-temp1*r12hat[2]);

  temp1 = -eta*u_r;
  temp2 = eta*chi;
  fforce[0] = temp1*dchi[0]-temp2*dUr[0];
  fforce[1] = temp1*dchi[1]-temp2*dUr[1];
  fforce[2] = temp1*dchi[2]-temp2*dUr[2];
  temp3 = eta*chi;
  fforce[0] = temp1*dchi[0]-temp3*dUr[0];
  fforce[1] = temp1*dchi[1]-temp3*dUr[1];
  fforce[2] = temp1*dchi[2]-temp3*dUr[2];

  // torque for particle 1 and 2
  // compute dUr
@@ -668,18 +668,17 @@ double PairGayBerne::gayberne_analytic(const int i,const int j,double a1[3][3],

  MathExtra::vecmat(iota,b1,tempv);
  MathExtra::cross3(tempv,iota,dchi);
  temp1 = -4.0/rsq;
  dchi[0] *= temp1;
  dchi[1] *= temp1;
  dchi[2] *= temp1;
  dchi[0] *= temp2;
  dchi[1] *= temp2;
  dchi[2] *= temp2;
  double dchi2[3];

  if (newton_pair || j < nlocal) {
    MathExtra::vecmat(iota,b2,tempv);
    MathExtra::cross3(tempv,iota,dchi2);
    dchi2[0] *= temp1;
    dchi2[1] *= temp1;
    dchi2[2] *= temp1;
    dchi2[0] *= temp2;
    dchi2[1] *= temp2;
    dchi2[2] *= temp2;
  }

  // compute d_eta
+10 −11
Original line number Diff line number Diff line
@@ -555,10 +555,10 @@ void PairGayBerneIntel::eval(const int offload, const int vflag,
          dchi_2 = temp2 * (iota_2 - temp1 * r12hat_2);

          temp1 = -eta * u_r;
          temp2 = eta * chi;
          fforce_0 = temp1 * dchi_0 - temp2 * dUr_0;
          fforce_1 = temp1 * dchi_1 - temp2 * dUr_1;
          fforce_2 = temp1 * dchi_2 - temp2 * dUr_2;
          temp3 = eta * chi;
          fforce_0 = temp1 * dchi_0 - temp3 * dUr_0;
          fforce_1 = temp1 * dchi_1 - temp3 * dUr_1;
          fforce_2 = temp1 * dchi_2 - temp3 * dUr_2;

          // torque for particle 1 and 2
          // compute dUr
@@ -579,18 +579,17 @@ void PairGayBerneIntel::eval(const int offload, const int vflag,

          ME_vecmat(iota, b1, tempv);
          ME_cross3(tempv, iota, dchi);
          temp1 = (flt_t)-4.0 / rsq_form[jj];
          dchi_0 *= temp1;
          dchi_1 *= temp1;
          dchi_2 *= temp1;
          dchi_0 *= temp2;
          dchi_1 *= temp2;
          dchi_2 *= temp2;
          flt_t dchi2_0, dchi2_1, dchi2_2;

          if (NEWTON_PAIR) {
            ME_vecmat(iota, b2, tempv);
            ME_cross3(tempv, iota, dchi2);
            dchi2_0 *= temp1;
            dchi2_1 *= temp1;
            dchi2_2 *= temp1;
            dchi2_0 *= temp2;
            dchi2_1 *= temp2;
            dchi2_2 *= temp2;
          }

          // compute d_eta