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 537 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.

00538 {
00539   enum { D = CMAP_TABLE_DIM };
00540   enum { N = CMAP_SETUP_DIM };
00541   const double h_1 = 1.0 / ( CMAP_SPACING * PI / 180.0) ;
00542   const double h_2 = h_1 * h_1;
00543   const double tr_h = 3.0 * h_1;
00544   int i, j;
00545   int ij;
00546   int ijp1p1, ijm1m1, ijp1m1, ijm1p1;
00547   int ijp2p2, ijm2m2, ijp2m2, ijm2p2;
00548 
00549   /* allocate spline coefficient matrix */
00550   double* const m = new double[N*N];
00551   memset(m,0,N*N*sizeof(double));
00552 
00553   /* initialize spline coefficient matrix */
00554   m[0] = 4;
00555   for (i = 1;  i < N;  i++) {
00556     m[INDEX(N,i-1,i)] = 1;
00557     m[INDEX(N,i,i-1)] = 1;
00558     m[INDEX(N,i,i)] = 4;
00559   }
00560   /* periodic boundary conditions for spline */
00561   m[INDEX(N,0,N-1)] = 1;
00562   m[INDEX(N,N-1,0)] = 1;
00563 
00564   /* compute LU-decomposition for this matrix */
00565   lu_decomp_nopivot(m, N);
00566 
00567   /* allocate vector for solving spline derivatives */
00568   double* const v = new double[N];
00569   memset(v,0,N*sizeof(double));
00570 
00571     /* march through rows of table */
00572     for (i = 0;  i < N;  i++) {
00573 
00574       /* setup RHS vector for solving spline derivatives */
00575       v[0] = tr_h * (table[INDEX(D,i,1)].d00 - table[INDEX(D,i,N-1)].d00);
00576       for (j = 1;  j < N;  j++) {
00577         v[j] = tr_h * (table[INDEX(D,i,j+1)].d00 - table[INDEX(D,i,j-1)].d00);
00578       }
00579 
00580       /* solve system, returned into vector */
00581       forward_back_sub(v, m, N);
00582 
00583       /* store values as derivatives wrt differenced table values */
00584       for (j = 0;  j < N;  j++) {
00585         table[INDEX(D,i,j)].d01 = v[j];
00586       }
00587       table[INDEX(D,i,N)].d01 = v[0];
00588     }
00589     for (j = 0;  j <= N;  j++) {
00590       table[INDEX(D,N,j)].d01 = table[INDEX(D,0,j)].d01;
00591     }
00592 
00593     /* march through columns of table */
00594     for (j = 0;  j < N;  j++) {
00595 
00596       /* setup RHS vector for solving spline derivatives */
00597       v[0] = tr_h * (table[INDEX(D,1,j)].d00 - table[INDEX(D,N-1,j)].d00);
00598       for (i = 1;  i < N;  i++) {
00599         v[i] = tr_h * (table[INDEX(D,i+1,j)].d00 - table[INDEX(D,i-1,j)].d00);
00600       }
00601 
00602       /* solve system, returned into vector */
00603       forward_back_sub(v, m, N);
00604 
00605       /* store values as derivatives wrt differenced table values */
00606       for (i = 0;  i < N;  i++) {
00607         table[INDEX(D,i,j)].d10 = v[i];
00608       }
00609       table[INDEX(D,N,j)].d10 = v[0];
00610     }
00611     for (i = 0;  i <= N;  i++) {
00612       table[INDEX(D,i,N)].d10 = table[INDEX(D,i,0)].d10;
00613     }
00614 
00615     /* march back through rows of table
00616      *
00617      * This is CHARMM's approach for calculating mixed partial derivatives,
00618      * by splining the first derivative values.
00619      *
00620      * Here we spline the dx values along y to calculate dxdy derivatives.
00621      *
00622      * Test cases show error with CHARMM is within 1e-5 roundoff error.
00623      */
00624     for (i = 0;  i < N;  i++) {
00625 
00626       /* setup RHS vector for solving dxy derivatives from dx */
00627       v[0] = tr_h * (table[INDEX(D,i,1)].d10 - table[INDEX(D,i,N-1)].d10);
00628       for (j = 1;  j < N;  j++) {
00629         v[j] = tr_h * (table[INDEX(D,i,j+1)].d10 - table[INDEX(D,i,j-1)].d10);
00630       }
00631 
00632       /* solve system, returned into vector */
00633       forward_back_sub(v, m, N);
00634 
00635       /* store values as dxy derivatives wrt differenced table values */
00636       for (j = 0;  j < N;  j++) {
00637         table[INDEX(D,i,j)].d11 = v[j];
00638       }
00639       table[INDEX(D,i,N)].d11 = v[0];
00640     }
00641     for (j = 0;  j <= N;  j++) {
00642       table[INDEX(D,N,j)].d11 = table[INDEX(D,0,j)].d11;
00643     }
00644 
00645   /* done with temp storage */
00646   delete [] m;
00647   delete [] v;
00648 
00649 }

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

Definition at line 669 of file ComputeCrossterms.C.

References INDEX, and j.

Referenced by crossterm_setup().

00670 {
00671   int i, j;
00672 
00673   /* in place forward elimination (solves Ly=b) using lower triangle of m */
00674   for (j = 0;  j < n-1;  j++) {
00675     for (i = j+1;  i < n;  i++) {
00676       b[i] -= m[INDEX(n,i,j)] * b[j];
00677     }
00678   }
00679   /* in place back substitution (solves Ux=y) using upper triangle of m */
00680   for (j = n-1;  j >= 0;  j--) {
00681     b[j] /= m[INDEX(n,j,j)];
00682     for (i = j-1;  i >= 0;  i--) {
00683       b[i] -= m[INDEX(n,i,j)] * b[j];
00684     }
00685   }
00686 }

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

Definition at line 652 of file ComputeCrossterms.C.

References INDEX, and j.

Referenced by crossterm_setup().

00653 {
00654   double l_ik;
00655   int i, j, k;
00656 
00657   for (k = 0;  k < n-1;  k++) {
00658     for (i = k+1;  i < n;  i++) {
00659       l_ik = m[INDEX(n,i,k)] / m[INDEX(n,k,k)];
00660       for (j = k;  j < n;  j++) {
00661         m[INDEX(n,i,j)] -= l_ik * m[INDEX(n,k,j)];
00662       }
00663       m[INDEX(n,i,k)] = l_ik;
00664     }
00665   }
00666 }


Generated on Sat Sep 23 01:17:16 2017 for NAMD by  doxygen 1.4.7