ComputeCrossterms.C File Reference

#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"

Go to the source code of this file.

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

#define FASTER

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

Definition at line 17 of file ComputeCrossterms.C.

#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

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

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;

}

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];
    }
  }
}

void lu_decomp_nopivot ( double *  m, int  n  )

static

Definition at line 704 of file ComputeCrossterms.C.

References INDEX.

Referenced by crossterm_setup().

{
  double l_ik;
  int i, j, k;

  for (k = 0;  k < n-1;  k++) {
    for (i = k+1;  i < n;  i++) {
      l_ik = m[INDEX(n,i,k)] / m[INDEX(n,k,k)];
      for (j = k;  j < n;  j++) {
        m[INDEX(n,i,j)] -= l_ik * m[INDEX(n,k,j)];
      }
      m[INDEX(n,i,k)] = l_ik;
    }
  }
}