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.

Defines

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


Define 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(), lu_decomp_nopivot(), Parameters::read_energy_type_bothcubspline(), and Parameters::read_energy_type_cubspline().


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.

00019      {
00020   CMAP_TABLE_DIM = 25,
00021   CMAP_SETUP_DIM = 24,
00022   CMAP_SPACING = 15,
00023   CMAP_PHI_0 = -180,
00024   CMAP_PSI_0 = -180
00025 };


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::d01, CrosstermData::d10, CrosstermData::d11, forward_back_sub(), INDEX, j, lu_decomp_nopivot(), and PI.

00590 {
00591   enum { D = CMAP_TABLE_DIM };
00592   enum { N = CMAP_SETUP_DIM };
00593   const double h_1 = 1.0 / ( CMAP_SPACING * PI / 180.0) ;
00594   const double h_2 = h_1 * h_1;
00595   const double tr_h = 3.0 * h_1;
00596   int i, j;
00597   int ij;
00598   int ijp1p1, ijm1m1, ijp1m1, ijm1p1;
00599   int ijp2p2, ijm2m2, ijp2m2, ijm2p2;
00600 
00601   /* allocate spline coefficient matrix */
00602   double* const m = new double[N*N];
00603   memset(m,0,N*N*sizeof(double));
00604 
00605   /* initialize spline coefficient matrix */
00606   m[0] = 4;
00607   for (i = 1;  i < N;  i++) {
00608     m[INDEX(N,i-1,i)] = 1;
00609     m[INDEX(N,i,i-1)] = 1;
00610     m[INDEX(N,i,i)] = 4;
00611   }
00612   /* periodic boundary conditions for spline */
00613   m[INDEX(N,0,N-1)] = 1;
00614   m[INDEX(N,N-1,0)] = 1;
00615 
00616   /* compute LU-decomposition for this matrix */
00617   lu_decomp_nopivot(m, N);
00618 
00619   /* allocate vector for solving spline derivatives */
00620   double* const v = new double[N];
00621   memset(v,0,N*sizeof(double));
00622 
00623     /* march through rows of table */
00624     for (i = 0;  i < N;  i++) {
00625 
00626       /* setup RHS vector for solving spline derivatives */
00627       v[0] = tr_h * (table[INDEX(D,i,1)].d00 - table[INDEX(D,i,N-1)].d00);
00628       for (j = 1;  j < N;  j++) {
00629         v[j] = tr_h * (table[INDEX(D,i,j+1)].d00 - table[INDEX(D,i,j-1)].d00);
00630       }
00631 
00632       /* solve system, returned into vector */
00633       forward_back_sub(v, m, N);
00634 
00635       /* store values as derivatives wrt differenced table values */
00636       for (j = 0;  j < N;  j++) {
00637         table[INDEX(D,i,j)].d01 = v[j];
00638       }
00639       table[INDEX(D,i,N)].d01 = v[0];
00640     }
00641     for (j = 0;  j <= N;  j++) {
00642       table[INDEX(D,N,j)].d01 = table[INDEX(D,0,j)].d01;
00643     }
00644 
00645     /* march through columns of table */
00646     for (j = 0;  j < N;  j++) {
00647 
00648       /* setup RHS vector for solving spline derivatives */
00649       v[0] = tr_h * (table[INDEX(D,1,j)].d00 - table[INDEX(D,N-1,j)].d00);
00650       for (i = 1;  i < N;  i++) {
00651         v[i] = tr_h * (table[INDEX(D,i+1,j)].d00 - table[INDEX(D,i-1,j)].d00);
00652       }
00653 
00654       /* solve system, returned into vector */
00655       forward_back_sub(v, m, N);
00656 
00657       /* store values as derivatives wrt differenced table values */
00658       for (i = 0;  i < N;  i++) {
00659         table[INDEX(D,i,j)].d10 = v[i];
00660       }
00661       table[INDEX(D,N,j)].d10 = v[0];
00662     }
00663     for (i = 0;  i <= N;  i++) {
00664       table[INDEX(D,i,N)].d10 = table[INDEX(D,i,0)].d10;
00665     }
00666 
00667     /* march back through rows of table
00668      *
00669      * This is CHARMM's approach for calculating mixed partial derivatives,
00670      * by splining the first derivative values.
00671      *
00672      * Here we spline the dx values along y to calculate dxdy derivatives.
00673      *
00674      * Test cases show error with CHARMM is within 1e-5 roundoff error.
00675      */
00676     for (i = 0;  i < N;  i++) {
00677 
00678       /* setup RHS vector for solving dxy derivatives from dx */
00679       v[0] = tr_h * (table[INDEX(D,i,1)].d10 - table[INDEX(D,i,N-1)].d10);
00680       for (j = 1;  j < N;  j++) {
00681         v[j] = tr_h * (table[INDEX(D,i,j+1)].d10 - table[INDEX(D,i,j-1)].d10);
00682       }
00683 
00684       /* solve system, returned into vector */
00685       forward_back_sub(v, m, N);
00686 
00687       /* store values as dxy derivatives wrt differenced table values */
00688       for (j = 0;  j < N;  j++) {
00689         table[INDEX(D,i,j)].d11 = v[j];
00690       }
00691       table[INDEX(D,i,N)].d11 = v[0];
00692     }
00693     for (j = 0;  j <= N;  j++) {
00694       table[INDEX(D,N,j)].d11 = table[INDEX(D,0,j)].d11;
00695     }
00696 
00697   /* done with temp storage */
00698   delete [] m;
00699   delete [] v;
00700 
00701 }

void forward_back_sub ( double *  b,
double *  m,
int  n 
) [static]

Definition at line 721 of file ComputeCrossterms.C.

References INDEX, and j.

Referenced by crossterm_setup().

00722 {
00723   int i, j;
00724 
00725   /* in place forward elimination (solves Ly=b) using lower triangle of m */
00726   for (j = 0;  j < n-1;  j++) {
00727     for (i = j+1;  i < n;  i++) {
00728       b[i] -= m[INDEX(n,i,j)] * b[j];
00729     }
00730   }
00731   /* in place back substitution (solves Ux=y) using upper triangle of m */
00732   for (j = n-1;  j >= 0;  j--) {
00733     b[j] /= m[INDEX(n,j,j)];
00734     for (i = j-1;  i >= 0;  i--) {
00735       b[i] -= m[INDEX(n,i,j)] * b[j];
00736     }
00737   }
00738 }

void lu_decomp_nopivot ( double *  m,
int  n 
) [static]

Definition at line 704 of file ComputeCrossterms.C.

References INDEX, and j.

Referenced by crossterm_setup().

00705 {
00706   double l_ik;
00707   int i, j, k;
00708 
00709   for (k = 0;  k < n-1;  k++) {
00710     for (i = k+1;  i < n;  i++) {
00711       l_ik = m[INDEX(n,i,k)] / m[INDEX(n,k,k)];
00712       for (j = k;  j < n;  j++) {
00713         m[INDEX(n,i,j)] -= l_ik * m[INDEX(n,k,j)];
00714       }
00715       m[INDEX(n,i,k)] = l_ik;
00716     }
00717   }
00718 }


Generated on Sat May 26 01:17:15 2018 for NAMD by  doxygen 1.4.7