Settle.C File Reference

#include "Settle.h"
#include <string.h>
#include <math.h>
#include "InfoStream.h"
#include <cmath>
#include <iostream>

Go to the source code of this file.

Functions
void	settle1init (BigReal pmO, BigReal pmH, BigReal hhdist, BigReal ohdist, BigReal &mO, BigReal &mH, BigReal &mOrmT, BigReal &mHrmT, BigReal &ra, BigReal &rb, BigReal &rc, BigReal &rra)
	initialize cached water properties More...
int	settle1 (const Vector *ref, Vector *pos, Vector *vel, BigReal invdt, BigReal mOrmT, BigReal mHrmT, BigReal ra, BigReal rb, BigReal rc, BigReal rra)
	optimized settle1 algorithm, reuses water properties as much as possible More...
template<int veclen>
void	settle1_SIMD (const Vector *ref, Vector *pos, BigReal mOrmT, BigReal mHrmT, BigReal ra, BigReal rb, BigReal rc, BigReal rra)
template<int veclen>
void	rattlePair (const RattleParam *rattleParam, const BigReal *refx, const BigReal *refy, const BigReal *refz, BigReal *posx, BigReal *posy, BigReal *posz, bool &consFailure)
static BigReal	det_3by3 (BigReal A[4][4])
static void	swap_row (BigReal A[4][4], BigReal b[4], int r1, int r2)
static void	solve_4by4 (BigReal lambda [4], BigReal A[4][4], BigReal sigma[4])
static void	solveMatrix (BigReal lambda [4], BigReal A[4][4], BigReal sigma[4], int icnt)
static void	solveFullInverse (BigReal A[4][4], BigReal S[4][4], int icnt)
void	MSHAKEIterate (const int icnt, const RattleParam *rattleParam, const BigReal *refx, const BigReal *refy, const BigReal *refz, BigReal *posx, BigReal *posy, BigReal *posz, const BigReal tol2, const int maxiter, bool &done, bool &consFailure)
void	LINCS (const int icnt, const RattleParam *rattleParam, const BigReal *refx, const BigReal *refy, const BigReal *refz, BigReal *posx, BigReal *posy, BigReal *posz, const BigReal tol2, const int maxiter, bool &done, bool &consFailure)
void	rattleN (const int icnt, const RattleParam *rattleParam, const BigReal *refx, const BigReal *refy, const BigReal *refz, BigReal *posx, BigReal *posy, BigReal *posz, const BigReal tol2, const int maxiter, bool &done, bool &consFailure)
template void	rattlePair< 1 > (const RattleParam *rattleParam, const BigReal *refx, const BigReal *refy, const BigReal *refz, BigReal *posx, BigReal *posy, BigReal *posz, bool &consFailure)
template void	settle1_SIMD< 2 > (const Vector *ref, Vector *pos, BigReal mOrmT, BigReal mHrmT, BigReal ra, BigReal rb, BigReal rc, BigReal rra)
template void	settle1_SIMD< 1 > (const Vector *ref, Vector *pos, BigReal mOrmT, BigReal mHrmT, BigReal ra, BigReal rb, BigReal rc, BigReal rra)
static int	settlev (const Vector *pos, BigReal ma, BigReal mb, Vector *vel, BigReal dt, Tensor *virial)
int	settle2 (BigReal mO, BigReal mH, const Vector *pos, Vector *vel, BigReal dt, Tensor *virial)
void	settle1_SOA (const double *__restrict ref_x, const double *__restrict ref_y, const double *__restrict ref_z, double *__restrict pos_x, double *__restrict pos_y, double *__restrict pos_z, int numWaters, BigReal mOrmT, BigReal mHrmT, BigReal ra, BigReal rb, BigReal rc, BigReal rra)

Function Documentation

det_3by3()

static BigReal det_3by3 ( BigReal  A[4][4] )

inlinestatic

Definition at line 711 of file Settle.C.

Referenced by solveMatrix().

{
    return
    A[0][0]*(A[1][1]*A[2][2]-A[1][2]*A[2][1])-
    A[0][1]*(A[1][0]*A[2][2]-A[1][2]*A[2][0])+
    A[0][2]*(A[1][0]*A[2][1]-A[1][1]*A[2][0]);
}

LINCS()

void LINCS	(	const int	icnt,
		const RattleParam *	rattleParam,
		const BigReal *	refx,
		const BigReal *	refy,
		const BigReal *	refz,
		BigReal *	posx,
		BigReal *	posy,
		BigReal *	posz,
		const BigReal	tol2,
		const int	maxiter,
		bool &	done,
		bool &	consFailure
	)

Definition at line 936 of file Settle.C.

References RattleParam::dsq, RattleParam::ia, RattleParam::ib, RattleParam::rma, RattleParam::rmb, and solveFullInverse().

Referenced by HomePatch::rattle1(), and HomePatch::rattle1_SOA().

{
    BigReal pabx[4];
    BigReal rabx[4]; 
    BigReal paby[4]; 
    BigReal raby[4]; 
    BigReal pabz[4]; 
    BigReal rabz[4];
    BigReal drab[4];

    //check each constraint
    consFailure = false;
    done = true;
    int iter = 0;
    #pragma unroll 
    for(int i = 0; i < 4; ++i)
    {
        int a = rattleParam[i].ia;
        int b = rattleParam[i].ib;
        pabx[i] = posx[a] - posx[b];
        paby[i] = posy[a] - posy[b];
        pabz[i] = posz[a] - posz[b];
        rabx[i] = refx[a] - refx[b];
        raby[i] = refy[a] - refy[b];
        rabz[i] = refz[a] - refz[b];
        BigReal pabsq = pabx[i]*pabx[i] + paby[i]*paby[i] + pabz[i]*pabz[i];
        BigReal rabsq = rattleParam[i].dsq;
        if(i < icnt)
            drab[i] = 1./sqrt(rabx[i]*rabx[i]+raby[i]*raby[i]+rabz[i]*rabz[i]);
        else
            drab[i] = 0;
        BigReal diffsq = pabsq - rabsq;
        if ( fabs(diffsq) > (rabsq * tol2)  && i < icnt)
            done = false;
    }
    if(!done)
    {
        BigReal S[4][4];
        BigReal A[4][4];
        //construct S
        #pragma unroll
        for(int i = 0; i < 16; ++i)
        {
            int idy = i >> 2;
            int idx = i - (idy << 2);
            S[idy][idx] = (rabx[idx]*rabx[idy]+raby[idx]*raby[idy]+rabz[idx]*rabz[idy])*rattleParam[idx].rma
                         *drab[idx]*drab[idy];
        }
        BigReal r_unc[4];
        #pragma unroll
        for(int i = 0; i < 4; ++i)
        {
            if(i < icnt)
                S[i][i] += rattleParam[i].rmb/rattleParam[i].rma*S[i][i];
            r_unc[i] = (pabx[i]*rabx[i]+paby[i]*raby[i]+pabz[i]*rabz[i])*drab[i] -sqrt(rattleParam[i].dsq);//- 1./drab[i];
        }
        BigReal lambda[4] = {0,0,0,0};
        #pragma unroll
        for(int i = 0; i < 4; ++i)
            #pragma unroll
            for(int j = 0; j < 4; ++j)
                A[i][j]=0;
        solveFullInverse(A, S, icnt);
        #pragma unroll
        for(int i = 0; i < 4; ++i)
            #pragma unroll
            for(int j = 0; j < 4; ++j)
                lambda[i] += A[i][j] * r_unc[j];
        #pragma unroll
        for(int i = 0; i < 4; ++i)
        {
            int a = rattleParam[i].ia;
            int b = rattleParam[i].ib;
            BigReal k   = lambda[i];
            BigReal rma = rattleParam[i].rma*k*drab[i];
            BigReal rmb = rattleParam[i].rmb*k*drab[i];
            posx[a] = posx[a] - rabx[i]*rma;
            posy[a] = posy[a] - raby[i]*rma;
            posz[a] = posz[a] - rabz[i]*rma;
            posx[b] = posx[b] + rabx[i]*rmb;
            posy[b] = posy[b] + raby[i]*rmb;
            posz[b] = posz[b] + rabz[i]*rmb;
        }
        for(iter = 1; iter < maxiter; ++iter)
        {
            done = true;
            #pragma unroll
            for(int i = 0; i < 4; ++i)
            {
                int a = rattleParam[i].ia;
                int b = rattleParam[i].ib;
                pabx[i] = posx[a] - posx[b];
                paby[i] = posy[a] - posy[b];
                pabz[i] = posz[a] - posz[b];
                BigReal pabsq = pabx[i]*pabx[i] + paby[i]*paby[i] + pabz[i]*pabz[i];
                BigReal rabsq = rattleParam[i].dsq;
                BigReal diffsq = pabsq - rabsq;
                if ( fabs(diffsq) > (rabsq * tol2) && i < icnt)
                    done = false;
                r_unc[i] = (pabx[i]*rabx[i]+paby[i]*raby[i]+pabz[i]*rabz[i])*drab[i] -sqrt(2.*rabsq-pabsq);
            }
            if(done)
                break;
            
            //solveMatrix(lambda, S, r_unc);
            lambda[0] = 0;
            lambda[1] = 0;
            lambda[2] = 0;
            lambda[3] = 0;
            #pragma unroll
            for(int i = 0; i < 4; ++i)
                #pragma unroll
                for(int j = 0; j < 4; ++j)
                    lambda[i] += A[i][j] * r_unc[j];

            #pragma unroll
            for(int i = 0; i < 4; ++i)
            {
                int a = rattleParam[i].ia;
                int b = rattleParam[i].ib;
                BigReal k   = lambda[i];
                BigReal rma = rattleParam[i].rma*k*drab[i];
                BigReal rmb = rattleParam[i].rmb*k*drab[i];
                posx[a] = posx[a] - rabx[i]*rma;
                posy[a] = posy[a] - raby[i]*rma;
                posz[a] = posz[a] - rabz[i]*rma;
                posx[b] = posx[b] + rabx[i]*rmb;
                posy[b] = posy[b] + raby[i]*rmb;
                posz[b] = posz[b] + rabz[i]*rmb;
            }
        }
    }

    if(iter >= maxiter)
        consFailure = true;
}

MSHAKEIterate()

void MSHAKEIterate	(	const int	icnt,
		const RattleParam *	rattleParam,
		const BigReal *	refx,
		const BigReal *	refy,
		const BigReal *	refz,
		BigReal *	posx,
		BigReal *	posy,
		BigReal *	posz,
		const BigReal	tol2,
		const int	maxiter,
		bool &	done,
		bool &	consFailure
	)

Definition at line 830 of file Settle.C.

References RattleParam::dsq, RattleParam::ia, RattleParam::ib, RattleParam::rma, RattleParam::rmb, and solveMatrix().

Referenced by HomePatch::rattle1(), and HomePatch::rattle1_SOA().

{
    BigReal sigma[4], lambda[4];
    BigReal A[4][4];

    BigReal pabx[4];
    BigReal rabx[4];
    BigReal paby[4];
    BigReal raby[4];
    BigReal pabz[4];
    BigReal rabz[4];
    register int loop;
    consFailure = false;
    done = true;
    #pragma unroll
    for(int i = 0; i < 4; ++i)
    {
        int a = rattleParam[i].ia;
        int b = rattleParam[i].ib;
        pabx[i] = posx[a] - posx[b];
        paby[i] = posy[a] - posy[b];
        pabz[i] = posz[a] - posz[b];
        rabx[i] = refx[a] - refx[b];
        raby[i] = refy[a] - refy[b];
        rabz[i] = refz[a] - refz[b];
    }
    #pragma unroll
    for(int i = 0; i < 4; ++i)
    {
        BigReal pabsq = pabx[i]*pabx[i] + paby[i]*paby[i] + pabz[i]*pabz[i];
        BigReal rabsq = rattleParam[i].dsq;
        BigReal diffsq = pabsq - rabsq;
        sigma[i] = diffsq;
        if ( fabs(diffsq) > (rabsq * tol2)  && i < icnt)
            done = false;
    }
    for(loop = 0; loop < maxiter; ++loop)
    {
        if(!done)
        {
            //construct A
            #pragma unroll
            for(int j = 0; j < 4; ++j)
            {
                BigReal rma = rattleParam[j].rma;
                #pragma unroll
                for(int i = 0; i < 4; ++i)
                {
                    A[j][i] = 2.*(pabx[j]*rabx[i]+paby[j]*raby[i]+pabz[j]*rabz[i])*rma;
                }
                BigReal rmb = rattleParam[j].rmb;
                A[j][j] += 2.*(pabx[j]*rabx[j]+paby[j]*raby[j]+pabz[j]*rabz[j])*rmb;
                lambda[j] = 0.;
            }
            lambda[0] = 0;
            lambda[1] = 0;
            lambda[2] = 0;
            lambda[3] = 0;
            solveMatrix(lambda, A, sigma, icnt);
            #pragma unroll
            for(int i = 0; i < 4; ++i)
            {
                int a = rattleParam[i].ia;
                int b = rattleParam[i].ib;
                BigReal rma = rattleParam[i].rma * lambda[i];
                BigReal rmb = rattleParam[i].rmb * lambda[i];

                posx[a] -= rma * rabx[i];
                posy[a] -= rma * raby[i];
                posz[a] -= rma * rabz[i];
                posx[b] += rmb * rabx[i];
                posy[b] += rmb * raby[i];
                posz[b] += rmb * rabz[i];
                
            }
        }
        else
            break;
        done = true;
        #pragma unroll
        for(int i = 0; i < 4; ++i)
        {
            int a = rattleParam[i].ia;
            int b = rattleParam[i].ib;
            pabx[i] = posx[a] - posx[b];
            paby[i] = posy[a] - posy[b];
            pabz[i] = posz[a] - posz[b];
            BigReal pabsq = pabx[i]*pabx[i] + paby[i]*paby[i] + pabz[i]*pabz[i];
            BigReal rabsq = rattleParam[i].dsq;
            BigReal diffsq = pabsq - rabsq;
            sigma[i] = diffsq;
            if ( fabs(diffsq) > (rabsq * tol2) && i < icnt)
                done = false;
        }
    }
    if(loop == maxiter)
    {
        consFailure = true;
    }
}

rattleN()

void rattleN	(	const int	icnt,
		const RattleParam *	rattleParam,
		const BigReal *	refx,
		const BigReal *	refy,
		const BigReal *	refz,
		BigReal *	posx,
		BigReal *	posy,
		BigReal *	posz,
		const BigReal	tol2,
		const int	maxiter,
		bool &	done,
		bool &	consFailure
	)

Definition at line 1359 of file Settle.C.

References RattleParam::dsq, RattleParam::ia, RattleParam::ib, RattleParam::rma, and RattleParam::rmb.

Referenced by HomePatch::rattle1(), and HomePatch::rattle1_SOA().

                                 {

  for (int iter = 0; iter < maxiter; ++iter ) {
    done = true;
    consFailure = false;
    for (int i = 0; i < icnt; ++i ) {
      int a = rattleParam[i].ia;
      int b = rattleParam[i].ib;
      BigReal pabx = posx[a] - posx[b];
      BigReal paby = posy[a] - posy[b];
      BigReal pabz = posz[a] - posz[b];
      BigReal pabsq = pabx*pabx + paby*paby + pabz*pabz;
      BigReal rabsq = rattleParam[i].dsq;
      BigReal diffsq = rabsq - pabsq;
      if ( fabs(diffsq) > (rabsq * tol2) ) {
        BigReal rabx = refx[a] - refx[b];
        BigReal raby = refy[a] - refy[b];
        BigReal rabz = refz[a] - refz[b];
        BigReal rpab = rabx*pabx + raby*paby + rabz*pabz;
        if ( rpab < ( rabsq * 1.0e-6 ) ) {
          done = false;
          consFailure = true;
          continue;
        }
        BigReal rma = rattleParam[i].rma;
        BigReal rmb = rattleParam[i].rmb;
        BigReal gab = diffsq / ( 2.0 * ( rma + rmb ) * rpab );
        BigReal dpx = rabx * gab;
        BigReal dpy = raby * gab;
        BigReal dpz = rabz * gab;
        posx[a] += rma * dpx;
        posy[a] += rma * dpy;
        posz[a] += rma * dpz;
        posx[b] -= rmb * dpx;
        posy[b] -= rmb * dpy;
        posz[b] -= rmb * dpz;
        done = false;
      }
    }
    if ( done ) break;
  }

}

rattlePair()

template<int veclen>

void rattlePair	(	const RattleParam *	rattleParam,
		const BigReal *	refx,
		const BigReal *	refy,
		const BigReal *	refz,
		BigReal *	posx,
		BigReal *	posy,
		BigReal *	posz,
		bool &	consFailure
	)

Definition at line 554 of file Settle.C.

References RattleParam::dsq, RattleParam::ia, RattleParam::ib, RattleParam::rma, and RattleParam::rmb.

                                                                  {

  int a = rattleParam[0].ia;
  int b = rattleParam[0].ib;
  BigReal pabx = posx[a] - posx[b];
  BigReal paby = posy[a] - posy[b];
  BigReal pabz = posz[a] - posz[b];
  BigReal pabsq = pabx*pabx + paby*paby + pabz*pabz;
  BigReal rabsq = rattleParam[0].dsq;
  BigReal diffsq = rabsq - pabsq;
  BigReal rabx = refx[a] - refx[b];
  BigReal raby = refy[a] - refy[b];
  BigReal rabz = refz[a] - refz[b];

  BigReal refsq = rabx*rabx + raby*raby + rabz*rabz;

  BigReal rpab = rabx*pabx + raby*paby + rabz*pabz;

  BigReal rma = rattleParam[0].rma;
  BigReal rmb = rattleParam[0].rmb;

  BigReal gab;
  BigReal sqrtarg = rpab*rpab + refsq*diffsq;
  if ( sqrtarg < 0. ) {
    consFailure = true;
    gab = 0.;
  } else {
    consFailure = false;
    gab = (-rpab + sqrt(sqrtarg))/(refsq*(rma + rmb));
  }

  BigReal dpx = rabx * gab;
  BigReal dpy = raby * gab;
  BigReal dpz = rabz * gab;
  posx[a] += rma * dpx;
  posy[a] += rma * dpy;
  posz[a] += rma * dpz;
  posx[b] -= rmb * dpx;
  posy[b] -= rmb * dpy;
  posz[b] -= rmb * dpz;

}

rattlePair< 1 >()

template void rattlePair< 1 >	(	const RattleParam *	rattleParam,
		const BigReal *	refx,
		const BigReal *	refy,
		const BigReal *	refz,
		BigReal *	posx,
		BigReal *	posy,
		BigReal *	posz,
		bool &	consFailure
	)

Referenced by HomePatch::rattle1(), and HomePatch::rattle1_SOA().

settle1()

int settle1	(	const Vector *	ref,
		Vector *	pos,
		Vector *	vel,
		BigReal	invdt,
		BigReal	mOrmT,
		BigReal	mHrmT,
		BigReal	ra,
		BigReal	rb,
		BigReal	rc,
		BigReal	rra
	)

optimized settle1 algorithm, reuses water properties as much as possible

Definition at line 63 of file Settle.C.

References Vector::unit(), Vector::x, Vector::y, and Vector::z.

Referenced by HomePatch::rattle1old().

                                                      {
#if defined(__SSE2__) && ! defined(NAMD_DISABLE_SSE)
  // SSE acceleration of some of the costly parts of settle using
  // the Intel C/C++ classes.  This implementation uses the SSE units
  // less efficiency than is potentially possible, but in order to do
  // better, the settle algorithm will have to be vectorized and operate
  // on multiple waters at a time.  Doing so could give us the ability to
  // do two (double precison) or four (single precision) waters at a time.
  // This code achieves a modest speedup without the need to reorganize
  // the NAMD structure.  Once we have water molecules sorted in a single
  // block we can do far better.

  // vectors in the plane of the original positions
  Vector b0, c0;

  __m128d REF0xy = _mm_loadu_pd((double *) &ref[0].x);  // ref0.y and ref0.x
  __m128d REF1xy = _mm_loadu_pd((double *) &ref[1].x);  // ref1.y and ref1.x

  __m128d B0xy = _mm_sub_pd(REF1xy, REF0xy);
  _mm_storeu_pd((double *) &b0.x, B0xy);
  b0.z = ref[1].z - ref[0].z;

  __m128d REF2xy = _mm_loadu_pd((double *) &ref[2].x);  // ref2.y and ref2.x

  __m128d C0xy = _mm_sub_pd(REF2xy, REF0xy);
  _mm_storeu_pd((double *) &c0.x, C0xy);
  c0.z = ref[2].z - ref[0].z;

  // new center of mass
  // Vector d0 = pos[0] * mOrmT + ((pos[1] + pos[2]) * mHrmT);
  __align(16) Vector a1;
  __align(16) Vector b1;
  __align(16) Vector c1;
  __align(16) Vector d0;

  __m128d POS1xy = _mm_loadu_pd((double *) &pos[1].x);
  __m128d POS2xy = _mm_loadu_pd((double *) &pos[2].x);
  __m128d PMHrmTxy = _mm_mul_pd(_mm_add_pd(POS1xy, POS2xy), _mm_set1_pd(mHrmT));

  __m128d POS0xy = _mm_loadu_pd((double *) &pos[0].x);
  __m128d PMOrmTxy = _mm_mul_pd(POS0xy, _mm_set1_pd(mOrmT));
  __m128d D0xy = _mm_add_pd(PMOrmTxy, PMHrmTxy);

  d0.z = pos[0].z * mOrmT + ((pos[1].z + pos[2].z) * mHrmT);
  a1.z = pos[0].z - d0.z;
  b1.z = pos[1].z - d0.z;
  c1.z = pos[2].z - d0.z;

  __m128d A1xy = _mm_sub_pd(POS0xy, D0xy);
  _mm_store_pd((double *) &a1.x, A1xy); // must be aligned

  __m128d B1xy = _mm_sub_pd(POS1xy, D0xy);
  _mm_store_pd((double *) &b1.x, B1xy); // must be aligned

  __m128d C1xy = _mm_sub_pd(POS2xy, D0xy);
  _mm_store_pd((double *) &c1.x, C1xy); // must be aligned

  _mm_store_pd((double *) &d0.x, D0xy); // must be aligned
  
  // Vectors describing transformation from original coordinate system to
  // the 'primed' coordinate system as in the diagram.  
  Vector n0 = cross(b0, c0);
  Vector n1 = cross(a1, n0); 
  Vector n2 = cross(n0, n1); 
#else
  // vectors in the plane of the original positions
  Vector b0 = ref[1]-ref[0];
  Vector c0 = ref[2]-ref[0];
  
  // new center of mass
  Vector d0 = pos[0]*mOrmT + ((pos[1] + pos[2])*mHrmT);
 
  Vector a1 = pos[0] - d0;
  Vector b1 = pos[1] - d0;
  Vector c1 = pos[2] - d0;
  
  // Vectors describing transformation from original coordinate system to
  // the 'primed' coordinate system as in the diagram.  
  Vector n0 = cross(b0, c0);
  Vector n1 = cross(a1, n0); 
  Vector n2 = cross(n0, n1); 
#endif

#if defined(__SSE2__) && ! defined(NAMD_DISABLE_SSE) && ! defined(MISSING_mm_cvtsd_f64)
  __m128d l1 = _mm_set_pd(n0.x, n0.y);
  l1 = _mm_mul_pd(l1, l1);
  // n0.x^2 + n0.y^2
  double l1xy0 = _mm_cvtsd_f64(_mm_add_sd(l1, _mm_shuffle_pd(l1, l1, 1)));

  __m128d l3 = _mm_set_pd(n1.y, n1.z);
  l3 = _mm_mul_pd(l3, l3);
  // n1.y^2 + n1.z^2
  double l3yz1 = _mm_cvtsd_f64(_mm_add_sd(l3, _mm_shuffle_pd(l3, l3, 1)));

  __m128d l2 = _mm_set_pd(n1.x, n0.z);
  // len(n1)^2 and len(n0)^2 
  __m128d ts01 = _mm_add_pd(_mm_set_pd(l3yz1, l1xy0), _mm_mul_pd(l2, l2));

  __m128d l4 = _mm_set_pd(n2.x, n2.y);
  l4 = _mm_mul_pd(l4, l4);
  // n2.x^2 + n2.y^2
  double l4xy2 = _mm_cvtsd_f64(_mm_add_sd(l4, _mm_shuffle_pd(l4, l4, 1)));
  double ts2 = l4xy2 + (n2.z * n2.z);              // len(n2)^2

  double  invlens[4];
  // since rsqrt_nr() doesn't work with current compiler
  // this is the next best option 
  static const __m128d fvecd1p0 = _mm_set1_pd(1.0);

  // 1/len(n1) and 1/len(n0)
  __m128d invlen12 = _mm_div_pd(fvecd1p0, _mm_sqrt_pd(ts01));

  // invlens[0]=1/len(n0), invlens[1]=1/len(n1)
  _mm_storeu_pd(invlens, invlen12);

  n0 = n0 * invlens[0];

  // shuffle the order of operations around from the normal algorithm so
  // that we can double pump sqrt() with n2 and cosphi at the same time
  // these components are usually computed down in the canonical water block
  BigReal A1Z = n0 * a1;
  BigReal sinphi = A1Z * rra;
  BigReal tmp = 1.0-sinphi*sinphi;

  __m128d n2cosphi = _mm_sqrt_pd(_mm_set_pd(tmp, ts2));
  // invlens[2] = 1/len(n2), invlens[3] = cosphi
  _mm_storeu_pd(invlens+2, n2cosphi);

  n1 = n1 * invlens[1];
  n2 = n2 * (1.0 / invlens[2]);
  BigReal cosphi = invlens[3];

  b0 = Vector(n1*b0, n2*b0, n0*b0); // note: b0.z is never referenced again
  c0 = Vector(n1*c0, n2*c0, n0*c0); // note: c0.z is never referenced again
 
  b1 = Vector(n1*b1, n2*b1, n0*b1);
  c1 = Vector(n1*c1, n2*c1, n0*c1);

  // now we can compute positions of canonical water 
  BigReal sinpsi = (b1.z - c1.z)/(2.0*rc*cosphi);
  tmp = 1.0-sinpsi*sinpsi;
  BigReal cospsi = sqrt(tmp);
#else
  n0 = n0.unit();
  n1 = n1.unit();
  n2 = n2.unit();

  b0 = Vector(n1*b0, n2*b0, n0*b0); // note: b0.z is never referenced again
  c0 = Vector(n1*c0, n2*c0, n0*c0); // note: c0.z is never referenced again
 
  BigReal A1Z = n0 * a1;
  b1 = Vector(n1*b1, n2*b1, n0*b1);
  c1 = Vector(n1*c1, n2*c1, n0*c1);

  // now we can compute positions of canonical water 
  BigReal sinphi = A1Z * rra;
  BigReal tmp = 1.0-sinphi*sinphi;
  BigReal cosphi = sqrt(tmp);
  BigReal sinpsi = (b1.z - c1.z)/(2.0*rc*cosphi);
  tmp = 1.0-sinpsi*sinpsi;
  BigReal cospsi = sqrt(tmp);
#endif

  BigReal rbphi = -rb*cosphi;
  BigReal tmp1 = rc*sinpsi*sinphi;
  BigReal tmp2 = rc*sinpsi*cosphi;
 
  Vector a2(0, ra*cosphi, ra*sinphi);
  Vector b2(-rc*cospsi, rbphi - tmp1, -rb*sinphi + tmp2);
  Vector c2( rc*cosphi, rbphi + tmp1, -rb*sinphi - tmp2);

  // there are no a0 terms because we've already subtracted the term off 
  // when we first defined b0 and c0.
  BigReal alpha = b2.x*(b0.x - c0.x) + b0.y*b2.y + c0.y*c2.y;
  BigReal beta  = b2.x*(c0.y - b0.y) + b0.x*b2.y + c0.x*c2.y;
  BigReal gama  = b0.x*b1.y - b1.x*b0.y + c0.x*c1.y - c1.x*c0.y;
 
  BigReal a2b2 = alpha*alpha + beta*beta;
  BigReal sintheta = (alpha*gama - beta*sqrt(a2b2 - gama*gama))/a2b2;
  BigReal costheta = sqrt(1.0 - sintheta*sintheta);
  
#if 0
  Vector a3( -a2.y*sintheta, 
              a2.y*costheta,
              a2.z);
  Vector b3(b2.x*costheta - b2.y*sintheta,
              b2.x*sintheta + b2.y*costheta,
              b2.z);
  Vector c3(c2.x*costheta - c2.y*sintheta,
              c2.x*sintheta + c2.y*costheta,
              c2.z);
  
#else
  Vector a3( -a2.y*sintheta, 
              a2.y*costheta,
              A1Z);
  Vector b3(b2.x*costheta - b2.y*sintheta,
              b2.x*sintheta + b2.y*costheta,
              b1.z);
  Vector c3(-b2.x*costheta - c2.y*sintheta,
            -b2.x*sintheta + c2.y*costheta,
              c1.z);

#endif

  // undo the transformation; generate new normal vectors from the transpose.
  Vector m1(n1.x, n2.x, n0.x);
  Vector m2(n1.y, n2.y, n0.y);
  Vector m0(n1.z, n2.z, n0.z);

  pos[0] = Vector(a3*m1, a3*m2, a3*m0) + d0;
  pos[1] = Vector(b3*m1, b3*m2, b3*m0) + d0;
  pos[2] = Vector(c3*m1, c3*m2, c3*m0) + d0;

  // dt can be negative during startup!
  if (invdt != 0) {
    vel[0] = (pos[0]-ref[0])*invdt;
    vel[1] = (pos[1]-ref[1])*invdt;
    vel[2] = (pos[2]-ref[2])*invdt;
  }

  return 0;
}

settle1_SIMD()

template<int veclen>

void settle1_SIMD	(	const Vector *	ref,
		Vector *	pos,
		BigReal	mOrmT,
		BigReal	mHrmT,
		BigReal	ra,
		BigReal	rb,
		BigReal	rc,
		BigReal	rra
	)

Definition at line 293 of file Settle.C.

References Vector::x, Vector::y, and Vector::z.

                                       {

  BigReal ref0xt[veclen];
  BigReal ref0yt[veclen];
  BigReal ref0zt[veclen];
  BigReal ref1xt[veclen];
  BigReal ref1yt[veclen];
  BigReal ref1zt[veclen];
  BigReal ref2xt[veclen];
  BigReal ref2yt[veclen];
  BigReal ref2zt[veclen];

  BigReal pos0xt[veclen];
  BigReal pos0yt[veclen];
  BigReal pos0zt[veclen];
  BigReal pos1xt[veclen];
  BigReal pos1yt[veclen];
  BigReal pos1zt[veclen];
  BigReal pos2xt[veclen];
  BigReal pos2yt[veclen];
  BigReal pos2zt[veclen];

  for (int i=0;i < veclen;i++) {
    ref0xt[i] = ref[i*3+0].x;
    ref0yt[i] = ref[i*3+0].y;
    ref0zt[i] = ref[i*3+0].z;
    ref1xt[i] = ref[i*3+1].x;
    ref1yt[i] = ref[i*3+1].y;
    ref1zt[i] = ref[i*3+1].z;
    ref2xt[i] = ref[i*3+2].x;
    ref2yt[i] = ref[i*3+2].y;
    ref2zt[i] = ref[i*3+2].z;

    pos0xt[i] = pos[i*3+0].x;
    pos0yt[i] = pos[i*3+0].y;
    pos0zt[i] = pos[i*3+0].z;
    pos1xt[i] = pos[i*3+1].x;
    pos1yt[i] = pos[i*3+1].y;
    pos1zt[i] = pos[i*3+1].z;
    pos2xt[i] = pos[i*3+2].x;
    pos2yt[i] = pos[i*3+2].y;
    pos2zt[i] = pos[i*3+2].z;
  }

#pragma omp simd
  for (int i=0;i < veclen;i++) {

    BigReal ref0x = ref0xt[i];
    BigReal ref0y = ref0yt[i];
    BigReal ref0z = ref0zt[i];
    BigReal ref1x = ref1xt[i];
    BigReal ref1y = ref1yt[i];
    BigReal ref1z = ref1zt[i];
    BigReal ref2x = ref2xt[i];
    BigReal ref2y = ref2yt[i];
    BigReal ref2z = ref2zt[i];

    BigReal pos0x = pos0xt[i];
    BigReal pos0y = pos0yt[i];
    BigReal pos0z = pos0zt[i];
    BigReal pos1x = pos1xt[i];
    BigReal pos1y = pos1yt[i];
    BigReal pos1z = pos1zt[i];
    BigReal pos2x = pos2xt[i];
    BigReal pos2y = pos2yt[i];
    BigReal pos2z = pos2zt[i];

    // vectors in the plane of the original positions
    BigReal b0x = ref1x - ref0x;
    BigReal b0y = ref1y - ref0y;
    BigReal b0z = ref1z - ref0z;

    BigReal c0x = ref2x - ref0x;
    BigReal c0y = ref2y - ref0y;
    BigReal c0z = ref2z - ref0z;
    
    // new center of mass
    BigReal d0x = pos0x*mOrmT + ((pos1x + pos2x)*mHrmT);
    BigReal d0y = pos0y*mOrmT + ((pos1y + pos2y)*mHrmT);
    BigReal d0z = pos0z*mOrmT + ((pos1z + pos2z)*mHrmT);
   
    BigReal a1x = pos0x - d0x;
    BigReal a1y = pos0y - d0y;
    BigReal a1z = pos0z - d0z;

    BigReal b1x = pos1x - d0x;
    BigReal b1y = pos1y - d0y;
    BigReal b1z = pos1z - d0z;

    BigReal c1x = pos2x - d0x;
    BigReal c1y = pos2y - d0y;
    BigReal c1z = pos2z - d0z;
    
    // Vectors describing transformation from original coordinate system to
    // the 'primed' coordinate system as in the diagram.
    // n0 = b0 x c0
    BigReal n0x = b0y*c0z-c0y*b0z;
    BigReal n0y = c0x*b0z-b0x*c0z;
    BigReal n0z = b0x*c0y-c0x*b0y;

    // n1 = a1 x n0
    BigReal n1x = a1y*n0z-n0y*a1z;
    BigReal n1y = n0x*a1z-a1x*n0z;
    BigReal n1z = a1x*n0y-n0x*a1y;

    // n2 = n0 x n1
    BigReal n2x = n0y*n1z-n1y*n0z;
    BigReal n2y = n1x*n0z-n0x*n1z;
    BigReal n2z = n0x*n1y-n1x*n0y;

    // Normalize n0
    BigReal n0inv = 1.0/sqrt(n0x*n0x + n0y*n0y + n0z*n0z);
    n0x *= n0inv;
    n0y *= n0inv;
    n0z *= n0inv;

    BigReal n1inv = 1.0/sqrt(n1x*n1x + n1y*n1y + n1z*n1z);
    n1x *= n1inv;
    n1y *= n1inv;
    n1z *= n1inv;

    BigReal n2inv = 1.0/sqrt(n2x*n2x + n2y*n2y + n2z*n2z);
    n2x *= n2inv;
    n2y *= n2inv;
    n2z *= n2inv;

    //b0 = Vector(n1*b0, n2*b0, n0*b0); // note: b0.z is never referenced again
    BigReal n1b0 = n1x*b0x + n1y*b0y + n1z*b0z;
    BigReal n2b0 = n2x*b0x + n2y*b0y + n2z*b0z;

    //c0 = Vector(n1*c0, n2*c0, n0*c0); // note: c0.z is never referenced again
    BigReal n1c0 = n1x*c0x + n1y*c0y + n1z*c0z;
    BigReal n2c0 = n2x*c0x + n2y*c0y + n2z*c0z;
   
    BigReal A1Z = n0x*a1x + n0y*a1y + n0z*a1z;
    
    //b1 = Vector(n1*b1, n2*b1, n0*b1);
    BigReal n1b1 = n1x*b1x + n1y*b1y + n1z*b1z;
    BigReal n2b1 = n2x*b1x + n2y*b1y + n2z*b1z;
    BigReal n0b1 = n0x*b1x + n0y*b1y + n0z*b1z;

    //c1 = Vector(n1*c1, n2*c1, n0*c1);
    BigReal n1c1 = n1x*c1x + n1y*c1y + n1z*c1z;
    BigReal n2c1 = n2x*c1x + n2y*c1y + n2z*c1z;
    BigReal n0c1 = n0x*c1x + n0y*c1y + n0z*c1z;

    // now we can compute positions of canonical water 
    BigReal sinphi = A1Z * rra;
    BigReal tmp = 1.0-sinphi*sinphi;
    BigReal cosphi = sqrt(tmp);
    BigReal sinpsi = (n0b1 - n0c1)/(2.0*rc*cosphi);
    tmp = 1.0-sinpsi*sinpsi;
    BigReal cospsi = sqrt(tmp);

    BigReal rbphi = -rb*cosphi;
    BigReal tmp1 = rc*sinpsi*sinphi;
    BigReal tmp2 = rc*sinpsi*cosphi;
   
    //Vector a2(0, ra*cosphi, ra*sinphi);
    BigReal a2y = ra*cosphi;

    //Vector b2(-rc*cospsi, rbphi - tmp1, -rb*sinphi + tmp2);
    BigReal b2x = -rc*cospsi;
    BigReal b2y = rbphi - tmp1;

    //Vector c2( rc*cosphi, rbphi + tmp1, -rb*sinphi - tmp2);
    BigReal c2y = rbphi + tmp1;

    // there are no a0 terms because we've already subtracted the term off 
    // when we first defined b0 and c0.
    BigReal alpha = b2x*(n1b0 - n1c0) + n2b0*b2y + n2c0*c2y;
    BigReal beta  = b2x*(n2c0 - n2b0) + n1b0*b2y + n1c0*c2y;
    BigReal gama  = n1b0*n2b1 - n1b1*n2b0 + n1c0*n2c1 - n1c1*n2c0;
   
    BigReal a2b2 = alpha*alpha + beta*beta;
    BigReal sintheta = (alpha*gama - beta*sqrt(a2b2 - gama*gama))/a2b2;
    BigReal costheta = sqrt(1.0 - sintheta*sintheta);
    
    //Vector a3( -a2y*sintheta, 
    //            a2y*costheta,
    //            A1Z);
    BigReal a3x = -a2y*sintheta;
    BigReal a3y = a2y*costheta;
    BigReal a3z = A1Z;

    // Vector b3(b2x*costheta - b2y*sintheta,
    //             b2x*sintheta + b2y*costheta,
    //             n0b1);
    BigReal b3x = b2x*costheta - b2y*sintheta;
    BigReal b3y = b2x*sintheta + b2y*costheta;
    BigReal b3z = n0b1;

    // Vector c3(-b2x*costheta - c2y*sintheta,
    //           -b2x*sintheta + c2y*costheta,
    //             n0c1);
    BigReal c3x = -b2x*costheta - c2y*sintheta;
    BigReal c3y = -b2x*sintheta + c2y*costheta;
    BigReal c3z = n0c1;

    // undo the transformation; generate new normal vectors from the transpose.
    // Vector m1(n1.x, n2.x, n0.x);
    BigReal m1x = n1x;
    BigReal m1y = n2x;
    BigReal m1z = n0x;

    // Vector m2(n1.y, n2.y, n0.y);
    BigReal m2x = n1y;
    BigReal m2y = n2y;
    BigReal m2z = n0y;

    // Vector m0(n1.z, n2.z, n0.z);
    BigReal m0x = n1z;
    BigReal m0y = n2z;
    BigReal m0z = n0z;

    //pos[i*3+0] = Vector(a3*m1, a3*m2, a3*m0) + d0;
    pos0x = a3x*m1x + a3y*m1y + a3z*m1z + d0x;
    pos0y = a3x*m2x + a3y*m2y + a3z*m2z + d0y;
    pos0z = a3x*m0x + a3y*m0y + a3z*m0z + d0z;

    // pos[i*3+1] = Vector(b3*m1, b3*m2, b3*m0) + d0;
    pos1x = b3x*m1x + b3y*m1y + b3z*m1z + d0x;
    pos1y = b3x*m2x + b3y*m2y + b3z*m2z + d0y;
    pos1z = b3x*m0x + b3y*m0y + b3z*m0z + d0z;

    // pos[i*3+2] = Vector(c3*m1, c3*m2, c3*m0) + d0;
    pos2x = c3x*m1x + c3y*m1y + c3z*m1z + d0x;
    pos2y = c3x*m2x + c3y*m2y + c3z*m2z + d0y;
    pos2z = c3x*m0x + c3y*m0y + c3z*m0z + d0z;

    pos0xt[i] = pos0x;
    pos0yt[i] = pos0y;
    pos0zt[i] = pos0z;
    pos1xt[i] = pos1x;
    pos1yt[i] = pos1y;
    pos1zt[i] = pos1z;
    pos2xt[i] = pos2x;
    pos2yt[i] = pos2y;
    pos2zt[i] = pos2z;
  }

  for (int i=0;i < veclen;i++) {
    pos[i*3+0].x = pos0xt[i];
    pos[i*3+0].y = pos0yt[i];
    pos[i*3+0].z = pos0zt[i];
    pos[i*3+1].x = pos1xt[i];
    pos[i*3+1].y = pos1yt[i];
    pos[i*3+1].z = pos1zt[i];
    pos[i*3+2].x = pos2xt[i];
    pos[i*3+2].y = pos2yt[i];
    pos[i*3+2].z = pos2zt[i];
  }

}

settle1_SIMD< 1 >()

template void settle1_SIMD< 1 >	(	const Vector *	ref,
		Vector *	pos,
		BigReal	mOrmT,
		BigReal	mHrmT,
		BigReal	ra,
		BigReal	rb,
		BigReal	rc,
		BigReal	rra
	)

Referenced by HomePatch::rattle1().

settle1_SIMD< 2 >()

template void settle1_SIMD< 2 >	(	const Vector *	ref,
		Vector *	pos,
		BigReal	mOrmT,
		BigReal	mHrmT,
		BigReal	ra,
		BigReal	rb,
		BigReal	rc,
		BigReal	rra
	)

Referenced by HomePatch::rattle1().

settle1_SOA()

void settle1_SOA	(	const double *__restrict	ref_x,
		const double *__restrict	ref_y,
		const double *__restrict	ref_z,
		double *__restrict	pos_x,
		double *__restrict	pos_y,
		double *__restrict	pos_z,
		int	numWaters,
		BigReal	mOrmT,
		BigReal	mHrmT,
		BigReal	ra,
		BigReal	rb,
		BigReal	rc,
		BigReal	rra
	)

Definition at line 1487 of file Settle.C.

Referenced by HomePatch::rattle1_SOA().

      {
  for (int i=0;  i < numWaters;  i++) {
    BigReal ref0x = ref_x[3*i];
    BigReal ref0y = ref_y[3*i];
    BigReal ref0z = ref_z[3*i];
    BigReal ref1x = ref_x[3*i+1];
    BigReal ref1y = ref_y[3*i+1];
    BigReal ref1z = ref_z[3*i+1];
    BigReal ref2x = ref_x[3*i+2];
    BigReal ref2y = ref_y[3*i+2];
    BigReal ref2z = ref_z[3*i+2];

    BigReal pos0x = pos_x[3*i];
    BigReal pos0y = pos_y[3*i];
    BigReal pos0z = pos_z[3*i];
    BigReal pos1x = pos_x[3*i+1];
    BigReal pos1y = pos_y[3*i+1];
    BigReal pos1z = pos_z[3*i+1];
    BigReal pos2x = pos_x[3*i+2];
    BigReal pos2y = pos_y[3*i+2];
    BigReal pos2z = pos_z[3*i+2];

    // vectors in the plane of the original positions
    BigReal b0x = ref1x - ref0x;
    BigReal b0y = ref1y - ref0y;
    BigReal b0z = ref1z - ref0z;

    BigReal c0x = ref2x - ref0x;
    BigReal c0y = ref2y - ref0y;
    BigReal c0z = ref2z - ref0z;
    
    // new center of mass
    BigReal d0x = pos0x*mOrmT + ((pos1x + pos2x)*mHrmT);
    BigReal d0y = pos0y*mOrmT + ((pos1y + pos2y)*mHrmT);
    BigReal d0z = pos0z*mOrmT + ((pos1z + pos2z)*mHrmT);
   
    BigReal a1x = pos0x - d0x;
    BigReal a1y = pos0y - d0y;
    BigReal a1z = pos0z - d0z;

    BigReal b1x = pos1x - d0x;
    BigReal b1y = pos1y - d0y;
    BigReal b1z = pos1z - d0z;

    BigReal c1x = pos2x - d0x;
    BigReal c1y = pos2y - d0y;
    BigReal c1z = pos2z - d0z;
    
    // Vectors describing transformation from original coordinate system to
    // the 'primed' coordinate system as in the diagram.
    // n0 = b0 x c0
    BigReal n0x = b0y*c0z-c0y*b0z;
    BigReal n0y = c0x*b0z-b0x*c0z;
    BigReal n0z = b0x*c0y-c0x*b0y;

    // n1 = a1 x n0
    BigReal n1x = a1y*n0z-n0y*a1z;
    BigReal n1y = n0x*a1z-a1x*n0z;
    BigReal n1z = a1x*n0y-n0x*a1y;

    // n2 = n0 x n1
    BigReal n2x = n0y*n1z-n1y*n0z;
    BigReal n2y = n1x*n0z-n0x*n1z;
    BigReal n2z = n0x*n1y-n1x*n0y;

    // Normalize n0
    BigReal n0inv = 1.0/sqrt(n0x*n0x + n0y*n0y + n0z*n0z);
    n0x *= n0inv;
    n0y *= n0inv;
    n0z *= n0inv;

    BigReal n1inv = 1.0/sqrt(n1x*n1x + n1y*n1y + n1z*n1z);
    n1x *= n1inv;
    n1y *= n1inv;
    n1z *= n1inv;

    BigReal n2inv = 1.0/sqrt(n2x*n2x + n2y*n2y + n2z*n2z);
    n2x *= n2inv;
    n2y *= n2inv;
    n2z *= n2inv;

    //b0 = Vector(n1*b0, n2*b0, n0*b0); // note: b0.z is never referenced again
    BigReal n1b0 = n1x*b0x + n1y*b0y + n1z*b0z;
    BigReal n2b0 = n2x*b0x + n2y*b0y + n2z*b0z;

    //c0 = Vector(n1*c0, n2*c0, n0*c0); // note: c0.z is never referenced again
    BigReal n1c0 = n1x*c0x + n1y*c0y + n1z*c0z;
    BigReal n2c0 = n2x*c0x + n2y*c0y + n2z*c0z;
   
    BigReal A1Z = n0x*a1x + n0y*a1y + n0z*a1z;
    
    //b1 = Vector(n1*b1, n2*b1, n0*b1);
    BigReal n1b1 = n1x*b1x + n1y*b1y + n1z*b1z;
    BigReal n2b1 = n2x*b1x + n2y*b1y + n2z*b1z;
    BigReal n0b1 = n0x*b1x + n0y*b1y + n0z*b1z;

    //c1 = Vector(n1*c1, n2*c1, n0*c1);
    BigReal n1c1 = n1x*c1x + n1y*c1y + n1z*c1z;
    BigReal n2c1 = n2x*c1x + n2y*c1y + n2z*c1z;
    BigReal n0c1 = n0x*c1x + n0y*c1y + n0z*c1z;

    // now we can compute positions of canonical water 
    BigReal sinphi = A1Z * rra;
    BigReal tmp = 1.0-sinphi*sinphi;
    BigReal cosphi = sqrt(tmp);
    BigReal sinpsi = (n0b1 - n0c1)/(2.0*rc*cosphi);
    tmp = 1.0-sinpsi*sinpsi;
    BigReal cospsi = sqrt(tmp);

    BigReal rbphi = -rb*cosphi;
    BigReal tmp1 = rc*sinpsi*sinphi;
    BigReal tmp2 = rc*sinpsi*cosphi;
   
    //Vector a2(0, ra*cosphi, ra*sinphi);
    BigReal a2y = ra*cosphi;

    //Vector b2(-rc*cospsi, rbphi - tmp1, -rb*sinphi + tmp2);
    BigReal b2x = -rc*cospsi;
    BigReal b2y = rbphi - tmp1;

    //Vector c2( rc*cosphi, rbphi + tmp1, -rb*sinphi - tmp2);
    BigReal c2y = rbphi + tmp1;

    // there are no a0 terms because we've already subtracted the term off 
    // when we first defined b0 and c0.
    BigReal alpha = b2x*(n1b0 - n1c0) + n2b0*b2y + n2c0*c2y;
    BigReal beta  = b2x*(n2c0 - n2b0) + n1b0*b2y + n1c0*c2y;
    BigReal gama  = n1b0*n2b1 - n1b1*n2b0 + n1c0*n2c1 - n1c1*n2c0;
   
    BigReal a2b2 = alpha*alpha + beta*beta;
    BigReal sintheta = (alpha*gama - beta*sqrt(a2b2 - gama*gama))/a2b2;
    BigReal costheta = sqrt(1.0 - sintheta*sintheta);
    
    //Vector a3( -a2y*sintheta, 
    //            a2y*costheta,
    //            A1Z);
    BigReal a3x = -a2y*sintheta;
    BigReal a3y = a2y*costheta;
    BigReal a3z = A1Z;

    // Vector b3(b2x*costheta - b2y*sintheta,
    //             b2x*sintheta + b2y*costheta,
    //             n0b1);
    BigReal b3x = b2x*costheta - b2y*sintheta;
    BigReal b3y = b2x*sintheta + b2y*costheta;
    BigReal b3z = n0b1;

    // Vector c3(-b2x*costheta - c2y*sintheta,
    //           -b2x*sintheta + c2y*costheta,
    //             n0c1);
    BigReal c3x = -b2x*costheta - c2y*sintheta;
    BigReal c3y = -b2x*sintheta + c2y*costheta;
    BigReal c3z = n0c1;

    // undo the transformation; generate new normal vectors from the transpose.
    // Vector m1(n1.x, n2.x, n0.x);
    BigReal m1x = n1x;
    BigReal m1y = n2x;
    BigReal m1z = n0x;

    // Vector m2(n1.y, n2.y, n0.y);
    BigReal m2x = n1y;
    BigReal m2y = n2y;
    BigReal m2z = n0y;

    // Vector m0(n1.z, n2.z, n0.z);
    BigReal m0x = n1z;
    BigReal m0y = n2z;
    BigReal m0z = n0z;

    //pos[i*3+0] = Vector(a3*m1, a3*m2, a3*m0) + d0;
    pos0x = a3x*m1x + a3y*m1y + a3z*m1z + d0x;
    pos0y = a3x*m2x + a3y*m2y + a3z*m2z + d0y;
    pos0z = a3x*m0x + a3y*m0y + a3z*m0z + d0z;

    // pos[i*3+1] = Vector(b3*m1, b3*m2, b3*m0) + d0;
    pos1x = b3x*m1x + b3y*m1y + b3z*m1z + d0x;
    pos1y = b3x*m2x + b3y*m2y + b3z*m2z + d0y;
    pos1z = b3x*m0x + b3y*m0y + b3z*m0z + d0z;

    // pos[i*3+2] = Vector(c3*m1, c3*m2, c3*m0) + d0;
    pos2x = c3x*m1x + c3y*m1y + c3z*m1z + d0x;
    pos2y = c3x*m2x + c3y*m2y + c3z*m2z + d0y;
    pos2z = c3x*m0x + c3y*m0y + c3z*m0z + d0z;

    pos_x[3*i] = pos0x;
    pos_y[3*i] = pos0y;
    pos_z[3*i] = pos0z;
    pos_x[3*i+1] = pos1x;
    pos_y[3*i+1] = pos1y;
    pos_z[3*i+1] = pos1z;
    pos_x[3*i+2] = pos2x;
    pos_y[3*i+2] = pos2y;
    pos_z[3*i+2] = pos2z;
  }

} // settle1_SOA()

settle1init()

void settle1init	(	BigReal	pmO,
		BigReal	pmH,
		BigReal	hhdist,
		BigReal	ohdist,
		BigReal &	mO,
		BigReal &	mH,
		BigReal &	mOrmT,
		BigReal &	mHrmT,
		BigReal &	ra,
		BigReal &	rb,
		BigReal &	rc,
		BigReal &	rra
	)

initialize cached water properties

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 46 of file Settle.C.

Referenced by HomePatch::buildRattleList(), HomePatch::buildRattleList_SOA(), and HomePatch::rattle1old().

                                                         {
    BigReal rmT = 1.0 / (pmO+pmH+pmH);
    mO = pmO;
    mH = pmH;
    mOrmT = pmO * rmT;
    mHrmT = pmH * rmT;
    BigReal t1 = 0.5*pmO/pmH;
    rc = 0.5*hhdist;
    ra = sqrt(ohdist*ohdist-rc*rc)/(1.0+t1);
    rb = t1*ra;
    rra = 1.0 / ra;
}

settle2()

int settle2	(	BigReal	mO,
		BigReal	mH,
		const Vector *	pos,
		Vector *	vel,
		BigReal	dt,
		Tensor *	virial
	)

Definition at line 1473 of file Settle.C.

References settlev().

Referenced by HomePatch::minimize_rattle2(), and HomePatch::rattle2().

                                                           {

  settlev(pos, mO, mH, vel, dt, virial);
  return 0;
}

settlev()

static int settlev ( const Vector *  pos, BigReal  ma, BigReal  mb, Vector *  vel, BigReal  dt, Tensor *  virial  )

static

Definition at line 1420 of file Settle.C.

References outer(), and Vector::unit().

Referenced by settle2().

                                                               {
  
  Vector rAB = pos[1]-pos[0];
  Vector rBC = pos[2]-pos[1];
  Vector rCA = pos[0]-pos[2];
 
  Vector AB = rAB.unit();
  Vector BC = rBC.unit();
  Vector CA = rCA.unit();
  
  BigReal cosA = -AB * CA;
  BigReal cosB = -BC * AB;
  BigReal cosC = -CA * BC;

  BigReal vab = (vel[1]-vel[0])*AB;
  BigReal vbc = (vel[2]-vel[1])*BC;
  BigReal vca = (vel[0]-vel[2])*CA;

  BigReal mab = ma+mb;
  
  BigReal d = (2*mab*mab + 2*ma*mb*cosA*cosB*cosC - 2*mb*mb*cosA*cosA
               - ma*mab*(cosB*cosB + cosC*cosC))*0.5/mb;

  BigReal tab = (vab*(2*mab - ma*cosC*cosC) +
                vbc*(mb*cosC*cosA - mab*cosB) +
                vca*(ma*cosB*cosC - 2*mb*cosA))*ma/d;
            
  BigReal tbc = (vbc*(mab*mab - mb*mb*cosA*cosA) +
                vca*ma*(mb*cosA*cosB - mab*cosC) +
                vab*ma*(mb*cosC*cosA - mab*cosB))/d;
  
  BigReal tca = (vca*(2*mab - ma*cosB*cosB) +
                vab*(ma*cosB*cosC - 2*mb*cosA) +
                vbc*(mb*cosA*cosB - mab*cosC))*ma/d;
 
  Vector ga = tab*AB - tca*CA;
  Vector gb = tbc*BC - tab*AB;
  Vector gc = tca*CA - tbc*BC;
#if 0
  if (virial) {
    *virial += 0.5*outer(tab, rAB)/dt;
    *virial += 0.5*outer(tbc, rBC)/dt;
    *virial += 0.5*outer(tca, rCA)/dt;
  }
#endif
  vel[0] += (0.5/ma)*ga;
  vel[1] += (0.5/mb)*gb;
  vel[2] += (0.5/mb)*gc;

  return 0;
}

solve_4by4()

static void solve_4by4 ( BigReal  lambda[4], BigReal  A[4][4], BigReal  sigma[4]  )

static

Definition at line 736 of file Settle.C.

References swap_row().

Referenced by solveMatrix().

{

    #pragma unroll
    for (int k = 0; k < 4; ++k)
    {
        #ifdef PIVOT
        int piv_row = k;
        BigReal Max = A[k][k];

        for (int row = k + 1; row < 4; ++row)
        {
            if ((tmp = fabs(A(row, k))) > Max)
            {
                piv_row = row;
                Max = tmp;
            }
        }
        if(k != piv_row)
            swap_row(A, sigma, k, piv_row);
        #endif
        for (int row = k + 1; row < 4; ++row)
        {
            BigReal tmp = A[row][k]/ A[k][k];
            for (int col = k+1; col < 4; col++)
                A[row][col] -= tmp * A[k][col];
            A[row][k]  = 0.;
            sigma[row]-= tmp * sigma[k];
        }
    }
    for (int row = 3; row >= 0; --row)
    {
        BigReal tmp = sigma[row];
        for (int j = 3; j > row; --j)
            tmp -= lambda[j] * A[row][j];
        lambda[row] = tmp / A[row][row];
    }
}

solveFullInverse()

static void solveFullInverse ( BigReal  A[4][4], BigReal  S[4][4], int  icnt  )

static

Definition at line 816 of file Settle.C.

References solveMatrix().

Referenced by LINCS().

{
    BigReal x[4];
    for(int i = 0; i < icnt; ++i)
    {
        BigReal b[4] = {0, 0, 0, 0};
        b[i] = 1.;
        solveMatrix(x, S, b, icnt);
        #pragma unroll
        for(int j = 0; j < 4; ++j) 
            A[j][i] = x[j];
    }
}

solveMatrix()

static void solveMatrix ( BigReal  lambda[4], BigReal  A[4][4], BigReal  sigma[4], int  icnt  )

static

Definition at line 775 of file Settle.C.

References det_3by3(), and solve_4by4().

Referenced by MSHAKEIterate(), and solveFullInverse().

{
    switch(icnt)
    {
        case 1:
        {
            lambda[0] = sigma[0]/A[0][0];
            break;
        }
        case 2:
        {
     
            BigReal det=1./(A[0][0]*A[1][1]-A[0][1]*A[1][0]);
            lambda[0]  = ( A[1][1]*sigma[0]-A[0][1]*sigma[1])*det;
            lambda[1]  = (-A[1][0]*sigma[0]+A[0][0]*sigma[1])*det;
            break;
        }
        case 3:
        {
            BigReal det = 1./det_3by3(A);
            lambda[0] = det*((A[1][1]*A[2][2]-A[1][2]*A[2][1])*sigma[0]-
                             (A[0][1]*A[2][2]-A[0][2]*A[2][1])*sigma[1]+
                             (A[0][1]*A[1][2]-A[0][2]*A[1][1])*sigma[2]);

            lambda[1] = det*((A[1][2]*A[2][0]-A[1][0]*A[2][2])*sigma[0]+
                             (A[0][0]*A[2][2]-A[0][2]*A[2][0])*sigma[1]-
                             (A[0][0]*A[1][2]-A[0][2]*A[1][0])*sigma[2]);

            lambda[2] = det*((A[1][0]*A[2][1]-A[1][1]*A[2][0])*sigma[0]-
                             (A[0][0]*A[2][1]-A[0][1]*A[2][0])*sigma[1]+
                             (A[0][0]*A[1][1]-A[0][1]*A[1][0])*sigma[2]);
            break;
        }
        case 4:
        {
            solve_4by4(lambda, A, sigma);
            break;
        }
    }
}

swap_row()

static void swap_row ( BigReal  A[4][4], BigReal  b[4], int  r1, int  r2  )

static

Definition at line 719 of file Settle.C.

Referenced by solve_4by4().

{
    #pragma unroll
    for (int i = 0; i < 4; i++)
    {
        BigReal* p1 = &A[r1][i];
        BigReal* p2 = &A[r2][i];
        BigReal tmp = *p1;
        *p1 = *p2;
        *p2 = tmp;
    }
    BigReal tmp;
    tmp   = b[r1];
    b[r1] = b[r2];
    b[r2] = tmp;
}