NAMD
Macros | Enumerations | Functions
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,
 
)    ((i)*ncols + (j))

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.

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().

590 {
591  enum { D = CMAP_TABLE_DIM };
592  enum { N = CMAP_SETUP_DIM };
593  const double h_1 = 1.0 / ( CMAP_SPACING * PI / 180.0) ;
594  const double h_2 = h_1 * h_1;
595  const double tr_h = 3.0 * h_1;
596  int i, j;
597  int ij;
598  int ijp1p1, ijm1m1, ijp1m1, ijm1p1;
599  int ijp2p2, ijm2m2, ijp2m2, ijm2p2;
600 
601  /* allocate spline coefficient matrix */
602  double* const m = new double[N*N];
603  memset(m,0,N*N*sizeof(double));
604 
605  /* initialize spline coefficient matrix */
606  m[0] = 4;
607  for (i = 1; i < N; i++) {
608  m[INDEX(N,i-1,i)] = 1;
609  m[INDEX(N,i,i-1)] = 1;
610  m[INDEX(N,i,i)] = 4;
611  }
612  /* periodic boundary conditions for spline */
613  m[INDEX(N,0,N-1)] = 1;
614  m[INDEX(N,N-1,0)] = 1;
615 
616  /* compute LU-decomposition for this matrix */
617  lu_decomp_nopivot(m, N);
618 
619  /* allocate vector for solving spline derivatives */
620  double* const v = new double[N];
621  memset(v,0,N*sizeof(double));
622 
623  /* march through rows of table */
624  for (i = 0; i < N; i++) {
625 
626  /* setup RHS vector for solving spline derivatives */
627  v[0] = tr_h * (table[INDEX(D,i,1)].d00 - table[INDEX(D,i,N-1)].d00);
628  for (j = 1; j < N; j++) {
629  v[j] = tr_h * (table[INDEX(D,i,j+1)].d00 - table[INDEX(D,i,j-1)].d00);
630  }
631 
632  /* solve system, returned into vector */
633  forward_back_sub(v, m, N);
634 
635  /* store values as derivatives wrt differenced table values */
636  for (j = 0; j < N; j++) {
637  table[INDEX(D,i,j)].d01 = v[j];
638  }
639  table[INDEX(D,i,N)].d01 = v[0];
640  }
641  for (j = 0; j <= N; j++) {
642  table[INDEX(D,N,j)].d01 = table[INDEX(D,0,j)].d01;
643  }
644 
645  /* march through columns of table */
646  for (j = 0; j < N; j++) {
647 
648  /* setup RHS vector for solving spline derivatives */
649  v[0] = tr_h * (table[INDEX(D,1,j)].d00 - table[INDEX(D,N-1,j)].d00);
650  for (i = 1; i < N; i++) {
651  v[i] = tr_h * (table[INDEX(D,i+1,j)].d00 - table[INDEX(D,i-1,j)].d00);
652  }
653 
654  /* solve system, returned into vector */
655  forward_back_sub(v, m, N);
656 
657  /* store values as derivatives wrt differenced table values */
658  for (i = 0; i < N; i++) {
659  table[INDEX(D,i,j)].d10 = v[i];
660  }
661  table[INDEX(D,N,j)].d10 = v[0];
662  }
663  for (i = 0; i <= N; i++) {
664  table[INDEX(D,i,N)].d10 = table[INDEX(D,i,0)].d10;
665  }
666 
667  /* march back through rows of table
668  *
669  * This is CHARMM's approach for calculating mixed partial derivatives,
670  * by splining the first derivative values.
671  *
672  * Here we spline the dx values along y to calculate dxdy derivatives.
673  *
674  * Test cases show error with CHARMM is within 1e-5 roundoff error.
675  */
676  for (i = 0; i < N; i++) {
677 
678  /* setup RHS vector for solving dxy derivatives from dx */
679  v[0] = tr_h * (table[INDEX(D,i,1)].d10 - table[INDEX(D,i,N-1)].d10);
680  for (j = 1; j < N; j++) {
681  v[j] = tr_h * (table[INDEX(D,i,j+1)].d10 - table[INDEX(D,i,j-1)].d10);
682  }
683 
684  /* solve system, returned into vector */
685  forward_back_sub(v, m, N);
686 
687  /* store values as dxy derivatives wrt differenced table values */
688  for (j = 0; j < N; j++) {
689  table[INDEX(D,i,j)].d11 = v[j];
690  }
691  table[INDEX(D,i,N)].d11 = v[0];
692  }
693  for (j = 0; j <= N; j++) {
694  table[INDEX(D,N,j)].d11 = table[INDEX(D,0,j)].d11;
695  }
696 
697  /* done with temp storage */
698  delete [] m;
699  delete [] v;
700 
701 }
static void forward_back_sub(double *b, double *m, int n)
#define INDEX(ncols, i, j)
static void lu_decomp_nopivot(double *m, int n)
BigReal d11
Definition: Parameters.h:119
#define PI
Definition: common.h:83
BigReal d00
Definition: Parameters.h:119
BigReal d01
Definition: Parameters.h:119
BigReal d10
Definition: Parameters.h:119
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().

722 {
723  int i, j;
724 
725  /* in place forward elimination (solves Ly=b) using lower triangle of m */
726  for (j = 0; j < n-1; j++) {
727  for (i = j+1; i < n; i++) {
728  b[i] -= m[INDEX(n,i,j)] * b[j];
729  }
730  }
731  /* in place back substitution (solves Ux=y) using upper triangle of m */
732  for (j = n-1; j >= 0; j--) {
733  b[j] /= m[INDEX(n,j,j)];
734  for (i = j-1; i >= 0; i--) {
735  b[i] -= m[INDEX(n,i,j)] * b[j];
736  }
737  }
738 }
#define INDEX(ncols, i, j)
void lu_decomp_nopivot ( double *  m,
int  n 
)
static

Definition at line 704 of file ComputeCrossterms.C.

References INDEX.

Referenced by crossterm_setup().

705 {
706  double l_ik;
707  int i, j, k;
708 
709  for (k = 0; k < n-1; k++) {
710  for (i = k+1; i < n; i++) {
711  l_ik = m[INDEX(n,i,k)] / m[INDEX(n,k,k)];
712  for (j = k; j < n; j++) {
713  m[INDEX(n,i,j)] -= l_ik * m[INDEX(n,k,j)];
714  }
715  m[INDEX(n,i,k)] = l_ik;
716  }
717  }
718 }
#define INDEX(ncols, i, j)