?? trilinterp.cc
字號:
/* trilinterp.cc */#include <math.h>#ifdef sgi #include "trilinterp.h"#else extern "C" { #include "trilinterp.h" }#endif/* linear interpolation from l (when a=0) to h (when a=1)*//* (equal to (a*h)+((1-a)*l) )*/#define LERP(a,l,h) ((l)+(((h)-(l))*(a)))extern int ElecMap;float trilinterp( CONST_FLOAT tcoord[MAX_ATOMS][SPACE], CONST_FLOAT charge[MAX_ATOMS], CONST_INT type[MAX_ATOMS], CONST_INT total_atoms, CONST_FLOAT map[MAX_GRID_PTS][MAX_GRID_PTS][MAX_GRID_PTS][MAX_MAPS], CONST_FLOAT inv_spacing, float elec[MAX_ATOMS], float emap[MAX_ATOMS], CONST_FLOAT xlo, /**/ CONST_FLOAT ylo, /* float lo[SPACE] ) SLOWER */ CONST_FLOAT zlo ) /**//*** float tcoord[MAX_ATOMS][SPACE]; temporary coordinates** float charge[MAX_ATOMS]; partial atomic charges** int type[MAX_ATOMS]; atom type of each atom** int total_atoms; number of atoms** float map[MAX_GRID_PTS][MAX_GRID_PTS][MAX_GRID_PTS][MAX_MAPS]; ** intermolecular interaction energies** float inv_spacing; = 1/(grid point spacing, in Angstroms) ** float elec[MAX_ATOMS]; electrostatic energies, atom by atom** float emap[MAX_ATOMS]; intermolecular energies** float xlo,ylo,zlo; minimum coordinates in x,y,z*//******************************************************************************//* Name: trilinterp *//* Function: Trilinear interpolation of interaction energies from map[] *//* using the coordinates in tcoord[]. *//* Copyright: (C) 1994, TSRI *//*----------------------------------------------------------------------------*//* Authors: Garrett M. Morris, TSRI, Accelerated C version 2.2 *//* David Goodsell, UCLA, Original FORTRAN version 1.0 *//* Date: 10/06/94 *//*----------------------------------------------------------------------------*//* Inputs: tcoord, charge, type, total_atoms, map, inv_spacing, lo *//* Returns: total energy *//* Globals: MAX_ATOMS, SPACE, MAX_ATOMS, MAX_GRID_PTS, MAX_MAPS. *//*----------------------------------------------------------------------------*//* Modification Record *//* Date Inits Comments *//* 05/05/91 GMM Translated into C. *//* 01/03/94 GMM Optimized code by examining 'cc -S trilinterp.c' output. *//* 10/06/94 GMM Optional 10% gain in speed, using nearest point, not *//* trilinear interpolation. Compile with -DMINPOINT flag. *//******************************************************************************/{ double emaptotal, electotal; double u, v, w; double p0u, p0v, p0w; double p1u, p1v, p1w; int AtomType; /* atom type */ int u0, v0, w0; int u1, v1, w1; register int i; /* i-th atom */#ifdef MINPOINT int x,y,z; /*MINPOINT*/#else double e, m; #endif emaptotal = electotal = 0.; for (i=0; i<total_atoms; i++) { AtomType = type[i]; u1 = (u0 = (int) (u = (tcoord[i][X]-xlo) * inv_spacing)) + 1; p1u = 1. - (p0u = u - (double) u0); v1 = (v0 = (int) (v = (tcoord[i][Y]-ylo) * inv_spacing)) + 1; p1v = 1. - (p0v = v - (double) v0); w1 = (w0 = (int) (w = (tcoord[i][Z]-zlo) * inv_spacing)) + 1; p1w = 1. - (p0w = w - (double) w0);#ifdef MINPOINT x = (p0u < p1u)? u0 : u1; /*MINPOINT*/ y = (p0v < p1v)? v0 : v1; /*MINPOINT*/ z = (p0w < p1w)? w0 : w1; /*MINPOINT*/ /*MINPOINT*/ electotal += (elec[i] = map[z][y][x][ElecMap] * charge[i]); /*MINPOINT*/ emaptotal += (emap[i] = map[z][y][x][AtomType]); /*MINPOINT*/#else e = m = 0.; e += p1u * p1v * p1w * map[ w0 ][ v0 ][ u0 ][ElecMap]; m += p1u * p1v * p1w * map[ w0 ][ v0 ][ u0 ][AtomType]; m += p0u * p1v * p1w * map[ w0 ][ v0 ][ u1 ][AtomType]; e += p0u * p1v * p1w * map[ w0 ][ v0 ][ u1 ][ElecMap]; e += p1u * p0v * p1w * map[ w0 ][ v1 ][ u0 ][ElecMap]; m += p1u * p0v * p1w * map[ w0 ][ v1 ][ u0 ][AtomType]; m += p0u * p0v * p1w * map[ w0 ][ v1 ][ u1 ][AtomType]; e += p0u * p0v * p1w * map[ w0 ][ v1 ][ u1 ][ElecMap]; e += p1u * p1v * p0w * map[ w1 ][ v0 ][ u0 ][ElecMap]; m += p1u * p1v * p0w * map[ w1 ][ v0 ][ u0 ][AtomType]; m += p0u * p1v * p0w * map[ w1 ][ v0 ][ u1 ][AtomType]; e += p0u * p1v * p0w * map[ w1 ][ v0 ][ u1 ][ElecMap]; e += p1u * p0v * p0w * map[ w1 ][ v1 ][ u0 ][ElecMap]; m += p1u * p0v * p0w * map[ w1 ][ v1 ][ u0 ][AtomType]; m += p0u * p0v * p0w * map[ w1 ][ v1 ][ u1 ][AtomType]; e += p0u * p0v * p0w * map[ w1 ][ v1 ][ u1 ][ElecMap]; electotal += (elec[i] = e * charge[i]); emaptotal += (emap[i] = m); #endif /* not MINPOINT */ }/*for 0 <= i < total_atoms*/ return( (float)(electotal + emaptotal) );}/* End of Function *//*----------------------------------------------------------------------------*//* quicktrilinterp.c */float quicktrilinterp( CONST_FLOAT tcoord[MAX_ATOMS][SPACE], CONST_FLOAT charge[MAX_ATOMS], CONST_INT type[MAX_ATOMS], CONST_INT total_atoms, CONST_FLOAT map[MAX_GRID_PTS][MAX_GRID_PTS][MAX_GRID_PTS][MAX_MAPS], CONST_FLOAT inv_spacing, CONST_FLOAT xlo, CONST_FLOAT ylo, CONST_FLOAT zlo ){ double etotal; double u, v, w; double p0u, p0v, p0w; double p1u, p1v, p1w; int AtomType; int u0, v0, w0; int u1, v1, w1; register int i; /* i-th atom */#ifdef MINPOINT int x,y,z;#else double e, m; #endif etotal = 0.; for (i=0; i<total_atoms; i++) { AtomType = type[i]; u1 = (u0 = (int) (u = (tcoord[i][X]-xlo) * inv_spacing)) + 1; p1u = 1. - (p0u = u - (double) u0); v1 = (v0 = (int) (v = (tcoord[i][Y]-ylo) * inv_spacing)) + 1; p1v = 1. - (p0v = v - (double) v0); w1 = (w0 = (int) (w = (tcoord[i][Z]-zlo) * inv_spacing)) + 1; p1w = 1. - (p0w = w - (double) w0);#ifdef MINPOINT x = (p0u < p1u)? u0 : u1; y = (p0v < p1v)? v0 : v1; z = (p0w < p1w)? w0 : w1; etotal += map[z][y][x][ElecMap] * charge[i] + map[z][y][x][AtomType]; #else e = m = 0.; e += p1u * p1v * p1w * map[ w0 ][ v0 ][ u0 ][ElecMap]; m += p1u * p1v * p1w * map[ w0 ][ v0 ][ u0 ][AtomType]; m += p0u * p1v * p1w * map[ w0 ][ v0 ][ u1 ][AtomType]; e += p0u * p1v * p1w * map[ w0 ][ v0 ][ u1 ][ElecMap]; e += p1u * p0v * p1w * map[ w0 ][ v1 ][ u0 ][ElecMap]; m += p1u * p0v * p1w * map[ w0 ][ v1 ][ u0 ][AtomType]; m += p1u * p1v * p0w * map[ w1 ][ v0 ][ u0 ][AtomType]; e += p1u * p1v * p0w * map[ w1 ][ v0 ][ u0 ][ElecMap]; e += p0u * p0v * p1w * map[ w0 ][ v1 ][ u1 ][ElecMap]; m += p0u * p0v * p1w * map[ w0 ][ v1 ][ u1 ][AtomType]; m += p1u * p0v * p0w * map[ w1 ][ v1 ][ u0 ][AtomType]; e += p1u * p0v * p0w * map[ w1 ][ v1 ][ u0 ][ElecMap]; e += p0u * p1v * p0w * map[ w1 ][ v0 ][ u1 ][ElecMap]; m += p0u * p1v * p0w * map[ w1 ][ v0 ][ u1 ][AtomType]; m += p0u * p0v * p0w * map[ w1 ][ v1 ][ u1 ][AtomType]; e += p0u * p0v * p0w * map[ w1 ][ v1 ][ u1 ][ElecMap]; etotal += e * charge[i] + m; #endif /* not MINPOINT */ }/*for 0 <= i < total_atoms*/ return( (float)etotal );}/* EOF *//*----------------------------------------------------------------------------*/float outsidetrilinterp(CONST_FLOAT tcoord[MAX_ATOMS][SPACE], CONST_FLOAT charge[MAX_ATOMS], CONST_INT type[MAX_ATOMS], CONST_INT total_atoms, CONST_FLOAT map[MAX_GRID_PTS][MAX_GRID_PTS][MAX_GRID_PTS][MAX_MAPS], CONST_FLOAT inv_spacing, // float elec[MAX_ATOMS], // float emap[MAX_ATOMS], CONST_FLOAT xlo, CONST_FLOAT ylo, CONST_FLOAT zlo, CONST_FLOAT xhi, CONST_FLOAT yhi, CONST_FLOAT zhi, CONST_FLOAT xcen, CONST_FLOAT ycen, CONST_FLOAT zcen ){ double etotal, epenalty; double x, y, z; double u, v, w; double p0u, p0v, p0w; double p1u, p1v, p1w; int AtomType; int u0, v0, w0; int u1, v1, w1; register int i; /* i-th atom */#ifdef MINPOINT int x,y,z;#else double e, m; #endif etotal = 0.; for (i=0; i<total_atoms; i++) { x = tcoord[i][X]; y = tcoord[i][Y]; z = tcoord[i][Z]; if (is_out_grid(x,y,z)) { x -= xcen; y -= ycen; z -= zcen; // sqhypotenuse(x,y,z) is the square of the distance // from grid's centre to atom epenalty = sqhypotenuse(x,y,z) * ENERGYPENALTY; // etotal += (elec[i] = epenalty) + (emap[i] = epenalty); etotal += (epenalty) + (epenalty); } else { AtomType = type[i]; u1 = (u0 = (int) (u = (tcoord[i][X]-xlo) * inv_spacing)) + 1; p1u = 1. - (p0u = u - (double) u0); v1 = (v0 = (int) (v = (tcoord[i][Y]-ylo) * inv_spacing)) + 1; p1v = 1. - (p0v = v - (double) v0); w1 = (w0 = (int) (w = (tcoord[i][Z]-zlo) * inv_spacing)) + 1; p1w = 1. - (p0w = w - (double) w0); #ifdef MINPOINT x = (p0u < p1u)? u0 : u1; y = (p0v < p1v)? v0 : v1; z = (p0w < p1w)? w0 : w1; // etotal += (elec[i] = map[z][y][x][ElecMap] * charge[i]) + // (emap[i] = map[z][y][x][AtomType]); etotal += (map[z][y][x][ElecMap] * charge[i]) + (map[i] = map[z][y][x][AtomType]); #else e = m = 0.; e += p1u * p1v * p1w * map[ w0 ][ v0 ][ u0 ][ElecMap]; m += p1u * p1v * p1w * map[ w0 ][ v0 ][ u0 ][AtomType]; m += p0u * p1v * p1w * map[ w0 ][ v0 ][ u1 ][AtomType]; e += p0u * p1v * p1w * map[ w0 ][ v0 ][ u1 ][ElecMap]; e += p1u * p0v * p1w * map[ w0 ][ v1 ][ u0 ][ElecMap]; m += p1u * p0v * p1w * map[ w0 ][ v1 ][ u0 ][AtomType]; m += p1u * p1v * p0w * map[ w1 ][ v0 ][ u0 ][AtomType]; e += p1u * p1v * p0w * map[ w1 ][ v0 ][ u0 ][ElecMap]; e += p0u * p0v * p1w * map[ w0 ][ v1 ][ u1 ][ElecMap]; m += p0u * p0v * p1w * map[ w0 ][ v1 ][ u1 ][AtomType]; m += p1u * p0v * p0w * map[ w1 ][ v1 ][ u0 ][AtomType]; e += p1u * p0v * p0w * map[ w1 ][ v1 ][ u0 ][ElecMap]; e += p0u * p1v * p0w * map[ w1 ][ v0 ][ u1 ][ElecMap]; m += p0u * p1v * p0w * map[ w1 ][ v0 ][ u1 ][AtomType]; m += p0u * p0v * p0w * map[ w1 ][ v1 ][ u1 ][AtomType]; e += p0u * p0v * p0w * map[ w1 ][ v1 ][ u1 ][ElecMap]; // etotal += (elec[i] = e * charge[i]) + (emap[i] = m); etotal += (e * charge[i]) + (m); #endif /* not MINPOINT */ } /* inside grid */ }/*for 0 <= i < total_atoms*/ return( (float)etotal );}/* EOF */?? 快捷鍵說明
復制代碼
Ctrl + C
搜索代碼
Ctrl + F
全屏模式
F11
切換主題
Ctrl + Shift + D
顯示快捷鍵
?
增大字號
Ctrl + =
減小字號
Ctrl + -