#include "InfoStream.h"
#include "ComputeCrossterms.h"
#include "Molecule.h"
#include "Parameters.h"
#include "Node.h"
#include "ReductionMgr.h"
#include "Lattice.h"
#include "PressureProfile.h"
#include "Debug.h"

Macros
#define	FASTER

#define	INDEX(ncols, i, j) ((i)*ncols + (j))

Enumerations
enum	{ CMAP_TABLE_DIM = 25, CMAP_SETUP_DIM = 24, CMAP_SPACING = 15, CMAP_PHI_0 = -180, CMAP_PSI_0 = -180 }

Functions
static void	lu_decomp_nopivot (double *m, int n)

static void	forward_back_sub (double b, double m, int n)

void	crossterm_setup (CrosstermData *table)

Macro Definition Documentation

◆ FASTER

#define FASTER

Definition at line 17 of file ComputeCrossterms.C.

◆ INDEX

#define INDEX	(	ncols,
		i,
		j
	)	((i)*ncols + (j))

Definition at line 27 of file ComputeCrossterms.C.

Referenced by CrosstermElem::computeForce(), crossterm_setup(), forward_back_sub(), and lu_decomp_nopivot().

Enumeration Type Documentation

◆ anonymous enum

anonymous enum

Enumerator
CMAP_TABLE_DIM
CMAP_SETUP_DIM
CMAP_SPACING
CMAP_PHI_0
CMAP_PSI_0

Definition at line 19 of file ComputeCrossterms.C.

      {
   CMAP_TABLE_DIM = 25,
   CMAP_SETUP_DIM = 24,
   CMAP_SPACING = 15,
   CMAP_PHI_0 = -180,
   CMAP_PSI_0 = -180
 };

Function Documentation

◆ crossterm_setup()

void crossterm_setup ( CrosstermData * table )

Definition at line 589 of file ComputeCrossterms.C.

References CMAP_SETUP_DIM, CMAP_SPACING, CMAP_TABLE_DIM, CrosstermData::d00, CrosstermData::d01, CrosstermData::d10, CrosstermData::d11, forward_back_sub(), INDEX, lu_decomp_nopivot(), and PI.

Referenced by Parameters::read_parm().

 {
   enum { D = CMAP_TABLE_DIM };
   enum { N = CMAP_SETUP_DIM };
   const double h_1 = 1.0 / ( CMAP_SPACING * PI / 180.0) ;
   const double h_2 = h_1 * h_1;
   const double tr_h = 3.0 * h_1;
   int i, j;
   int ij;
   int ijp1p1, ijm1m1, ijp1m1, ijm1p1;
   int ijp2p2, ijm2m2, ijp2m2, ijm2p2;
 
   /* allocate spline coefficient matrix */
   double* const m = new double[N*N];
   memset(m,0,N*N*sizeof(double));
 
   /* initialize spline coefficient matrix */
   m[0] = 4;
   for (i = 1;  i < N;  i++) {
     m[INDEX(N,i-1,i)] = 1;
     m[INDEX(N,i,i-1)] = 1;
     m[INDEX(N,i,i)] = 4;
   }
   /* periodic boundary conditions for spline */
   m[INDEX(N,0,N-1)] = 1;
   m[INDEX(N,N-1,0)] = 1;
 
   /* compute LU-decomposition for this matrix */
   lu_decomp_nopivot(m, N);
 
   /* allocate vector for solving spline derivatives */
   double* const v = new double[N];
   memset(v,0,N*sizeof(double));
 
     /* march through rows of table */
     for (i = 0;  i < N;  i++) {
 
       /* setup RHS vector for solving spline derivatives */
       v[0] = tr_h * (table[INDEX(D,i,1)].d00 - table[INDEX(D,i,N-1)].d00);
       for (j = 1;  j < N;  j++) {
         v[j] = tr_h * (table[INDEX(D,i,j+1)].d00 - table[INDEX(D,i,j-1)].d00);
       }
 
       /* solve system, returned into vector */
       forward_back_sub(v, m, N);
 
       /* store values as derivatives wrt differenced table values */
       for (j = 0;  j < N;  j++) {
         table[INDEX(D,i,j)].d01 = v[j];
       }
       table[INDEX(D,i,N)].d01 = v[0];
     }
     for (j = 0;  j <= N;  j++) {
       table[INDEX(D,N,j)].d01 = table[INDEX(D,0,j)].d01;
     }
 
     /* march through columns of table */
     for (j = 0;  j < N;  j++) {
 
       /* setup RHS vector for solving spline derivatives */
       v[0] = tr_h * (table[INDEX(D,1,j)].d00 - table[INDEX(D,N-1,j)].d00);
       for (i = 1;  i < N;  i++) {
         v[i] = tr_h * (table[INDEX(D,i+1,j)].d00 - table[INDEX(D,i-1,j)].d00);
       }
 
       /* solve system, returned into vector */
       forward_back_sub(v, m, N);
 
       /* store values as derivatives wrt differenced table values */
       for (i = 0;  i < N;  i++) {
         table[INDEX(D,i,j)].d10 = v[i];
       }
       table[INDEX(D,N,j)].d10 = v[0];
     }
     for (i = 0;  i <= N;  i++) {
       table[INDEX(D,i,N)].d10 = table[INDEX(D,i,0)].d10;
     }
 
     /* march back through rows of table
      *
      * This is CHARMM's approach for calculating mixed partial derivatives,
      * by splining the first derivative values.
      *
      * Here we spline the dx values along y to calculate dxdy derivatives.
      *
      * Test cases show error with CHARMM is within 1e-5 roundoff error.
      */
     for (i = 0;  i < N;  i++) {
 
       /* setup RHS vector for solving dxy derivatives from dx */
       v[0] = tr_h * (table[INDEX(D,i,1)].d10 - table[INDEX(D,i,N-1)].d10);
       for (j = 1;  j < N;  j++) {
         v[j] = tr_h * (table[INDEX(D,i,j+1)].d10 - table[INDEX(D,i,j-1)].d10);
       }
 
       /* solve system, returned into vector */
       forward_back_sub(v, m, N);
 
       /* store values as dxy derivatives wrt differenced table values */
       for (j = 0;  j < N;  j++) {
         table[INDEX(D,i,j)].d11 = v[j];
       }
       table[INDEX(D,i,N)].d11 = v[0];
     }
     for (j = 0;  j <= N;  j++) {
       table[INDEX(D,N,j)].d11 = table[INDEX(D,0,j)].d11;
     }
 
   /* done with temp storage */
   delete [] m;
   delete [] v;
 
 }

◆ forward_back_sub()

void forward_back_sub	(	double *	b,
		double *	m,
		int	n
	)

static

Definition at line 721 of file ComputeCrossterms.C.

References INDEX.

Referenced by crossterm_setup().

 {
   int i, j;
 
   /* in place forward elimination (solves Ly=b) using lower triangle of m */
   for (j = 0;  j < n-1;  j++) {
     for (i = j+1;  i < n;  i++) {
       b[i] -= m[INDEX(n,i,j)] * b[j];
     }
   }
   /* in place back substitution (solves Ux=y) using upper triangle of m */
   for (j = n-1;  j >= 0;  j--) {
     b[j] /= m[INDEX(n,j,j)];
     for (i = j-1;  i >= 0;  i--) {
       b[i] -= m[INDEX(n,i,j)] * b[j];
     }
   }
 }

◆ lu_decomp_nopivot()

void lu_decomp_nopivot	(	double *	m,
		int	n
	)