/*************************************************************************** forcefunchs.c - description ------------------- begin : Fri Mar 18 12:25:42 2005 copyright : (C) 2002 by Cavalli Andrea email : cavalli@bioc.unizh.ch **************************************************************************/ /*************************************************************************** * * * This program is free software; you can redistribute it and/or modify * * it under the terms of the GNU General Public License as published by * * the Free Software Foundation; either version 2 of the License, or * * (at your option) any later version. * * * ***************************************************************************/ #define COOR_DIM 4 void force_hssn_d(const int * first, const int * second, const int * boundarry, const int natoms, const double * coor, double * force, const double c2, const double c_on2, const int * type_code, const int type_size, const double * vdw_data){ int f; int fpos; int s; int spos; int b1,b2; double fracX; double fracY; double fracZ; double dx,dy,dz; double ix,iy,iz; double jx,jy,jz; double tx,ty,tz; double d2; double d6; double d2inv; double d6inv; double d12inv; double frac; double nb6; double nb12; double A; int tc1,tc2; int data_pos; const double c2inv = 1.0/c2; const double c6inv = c2inv*c2inv*c2inv; const double c12inv = c6inv*c6inv; const double c18inv = c12inv*c6inv; const double one = 1.0; int i,j; for(i = 0;ic_on2){ sh1 = c2-d2; sh2 = (c2+two*d2-three*c_on2); sw = sh1*sh2; sw1= four*(sh1-sh2); frac = nbvdw*sw*d2inv-(twelve_inv*A*d12inv)*sw1; frac = frac*c2diffinv*sh1; } else frac = nbvdw*d2inv; tx = frac*dx; ty = frac*dy; tz = frac*dz; fracX = fracX - tx; fracY = fracY - ty; fracZ = fracZ - tz; force[spos] = force[spos] + tx; force[spos+1] = force[spos+1] + ty; force[spos+2] = force[spos+2] + tz; } } force[fpos] = force[fpos] + fracX; force[fpos+1] = force[fpos+1] + fracY; force[fpos+2] = force[fpos+2] + fracZ; } } void force_hsswn_f(const int * first, const int * second, const int * boundarry, const int natoms, const float * coor, float * force, const float c2, const float c_on2, const int * type_code, const int type_size, const float * vdw_data){ int f; int fpos; int s; int spos; int b1,b2; float fracX; float fracY; float fracZ; float dx,dy,dz; float ix,iy,iz; float jx,jy,jz; float tx,ty,tz; float d2; float d6; float d2inv; float d6inv; float d12inv; float sh1; float sh2; float sw; float sw1; float frac; float nbvdw; float A; int tc1,tc2; int data_pos; const float c2diffinv = 1.0/((c2-c_on2)*(c2-c_on2)*(c2-c_on2)); const float one = 1.0; const float two = 2.0; const float three = 3.0; const float four = 4.0; const float six_inv = 1.0/6.0; const float twelve_inv = 1.0/12.0; int i,j; for(i = 0;ic_on2){ sh1 = c2-d2; sh2 = (c2+two*d2-three*c_on2); sw = sh1*sh2; sw1= four*(sh1-sh2); frac = nbvdw*sw*d2inv-(twelve_inv*A*d12inv)*sw1; frac = frac*c2diffinv*sh1; } else frac = nbvdw*d2inv; tx = frac*dx; ty = frac*dy; tz = frac*dz; fracX = fracX - tx; fracY = fracY - ty; fracZ = fracZ - tz; force[spos] = force[spos] + tx; force[spos+1] = force[spos+1] + ty; force[spos+2] = force[spos+2] + tz; } } force[fpos] = force[fpos] + fracX; force[fpos+1] = force[fpos+1] + fracY; force[fpos+2] = force[fpos+2] + fracZ; } } /*************************************************************************** IMPLEMENTATION OF NB ENERGY Elec Functional NONE Elec CUTOFF NONE VDW Functional LJ VDW CUTOFF SHIFT SOLVATION NONE Sol Functional SolvNone Geometry PLAIN Compute virial no Float type double CPU CPUPlain ***************************************************************************/ void energy_force_hssn_d (const int *first, const int *second, const int *boundary, const int natoms, const double *coor, double *force, const double c2, const double c_on2, const int *type_code, const int type_size, const double *vdw_data, double *ene_vdw, double *ene_cul) { int i; int j; int f; int fpos; int s; int spos; int b1; int b2; int tc1; int tc2; int data_pos; double fracX; double fracY; double fracZ; double dx, dy, dz; double ix, iy, iz; double jx, jy, jz; double tx, ty, tz; double d2; double d6; double d2inv; double d6inv; double d12inv; double vsh; double frac; double fvdw; double nb6; double nb12; double A; double evdw = 0; const double c2inv = 1.0 / c2; const double c6inv = c2inv * c2inv * c2inv; const double c12inv = c6inv * c6inv; const double c18inv = c12inv * c6inv; const double one = 1.0; const double six = 6.0; const double twelve = 12.0; /* Start outer loop over neighborlists */ for (i = 0; (i < natoms); i++) { /* Initialize outer atom */ f = first[i]; fpos = COOR_DIM * f; ix = coor[fpos]; iy = coor[fpos + 1]; iz = coor[fpos + 2]; /* Load limits for loop over neighbors */ b1 = boundary[2 * i]; b2 = boundary[2 * i + 1]; /* Zero forces */ fracX = 0; fracY = 0; fracZ = 0; /* Load pars */ tc1 = type_code[f]; /* Inner loop */ for (j = b1; (j < b2); j++) { /* Load coor second */ s = second[j]; spos = COOR_DIM * s; jx = coor[spos]; jy = coor[spos + 1]; jz = coor[spos + 2]; /* Compute distances */ dx = ix - jx; dy = iy - jy; dz = iz - jz; d2 = dx * dx + dy * dy + dz * dz; /* Compute force */ if (d2 < c2) { tc2 = type_code[s]; data_pos = 4 * (type_size * tc1 + tc2); A = vdw_data[data_pos]; d2inv = one / d2; d6 = d2 * d2 * d2; d6inv = d2inv * d2inv * d2inv; d12inv = d6inv * d6inv; vsh = one - d6 * c6inv; evdw = evdw + A * (d12inv - c12inv * (one + vsh + vsh)); nb12 = twelve * A * (d12inv - d6 * c18inv); fvdw = (nb12) * d2inv; frac = fvdw;; tx = frac * dx; ty = frac * dy; tz = frac * dz; fracX = fracX - tx; fracY = fracY - ty; fracZ = fracZ - tz; force[spos] = force[spos] + tx; force[spos + 1] = force[spos + 1] + ty; force[spos + 2] = force[spos + 2] + tz; } } /* Assign to first */ force[fpos] = force[fpos] + fracX; force[fpos + 1] = force[fpos + 1] + fracY; force[fpos + 2] = force[fpos + 2] + fracZ; } *ene_vdw = *ene_vdw + evdw; } /*************************************************************************** IMPLEMENTATION OF NB ENERGY Elec Functional NONE Elec CUTOFF NONE VDW Functional LJ VDW CUTOFF SHIFT SOLVATION NONE Sol Functional SolvNone Geometry PLAIN Compute virial no Float type float CPU CPUPlain ***************************************************************************/ void energy_force_hssn_f (const int *first, const int *second, const int *boundary, const int natoms, const float *coor, float *force, const float c2, const float c_on2, const int *type_code, const int type_size, const float *vdw_data, float *ene_vdw, float *ene_cul) { int i; int j; int f; int fpos; int s; int spos; int b1; int b2; int tc1; int tc2; int data_pos; float fracX; float fracY; float fracZ; float dx, dy, dz; float ix, iy, iz; float jx, jy, jz; float tx, ty, tz; float d2; float d6; float d2inv; float d6inv; float d12inv; float vsh; float frac; float fvdw; float nb6; float nb12; float A; float evdw = 0; const float c2inv = 1.0 / c2; const float c6inv = c2inv * c2inv * c2inv; const float c12inv = c6inv * c6inv; const float c18inv = c12inv * c6inv; const float one = 1.0; const float six = 6.0; const float twelve = 12.0; /* Start outer loop over neighborlists */ for (i = 0; (i < natoms); i++) { /* Initialize outer atom */ f = first[i]; fpos = COOR_DIM * f; ix = coor[fpos]; iy = coor[fpos + 1]; iz = coor[fpos + 2]; /* Load limits for loop over neighbors */ b1 = boundary[2 * i]; b2 = boundary[2 * i + 1]; /* Zero forces */ fracX = 0; fracY = 0; fracZ = 0; /* Load pars */ tc1 = type_code[f]; /* Inner loop */ for (j = b1; (j < b2); j++) { /* Load coor second */ s = second[j]; spos = COOR_DIM * s; jx = coor[spos]; jy = coor[spos + 1]; jz = coor[spos + 2]; /* Compute distances */ dx = ix - jx; dy = iy - jy; dz = iz - jz; d2 = dx * dx + dy * dy + dz * dz; /* Compute force */ if (d2 < c2) { tc2 = type_code[s]; data_pos = 4 * (type_size * tc1 + tc2); A = vdw_data[data_pos]; d2inv = one / d2; d6 = d2 * d2 * d2; d6inv = d2inv * d2inv * d2inv; d12inv = d6inv * d6inv; vsh = one - d6 * c6inv; evdw = evdw + A * (d12inv - c12inv * (one + vsh + vsh)); nb12 = twelve * A * (d12inv - d6 * c18inv); fvdw = (nb12) * d2inv; frac = fvdw;; tx = frac * dx; ty = frac * dy; tz = frac * dz; fracX = fracX - tx; fracY = fracY - ty; fracZ = fracZ - tz; force[spos] = force[spos] + tx; force[spos + 1] = force[spos + 1] + ty; force[spos + 2] = force[spos + 2] + tz; } } /* Assign to first */ force[fpos] = force[fpos] + fracX; force[fpos + 1] = force[fpos + 1] + fracY; force[fpos + 2] = force[fpos + 2] + fracZ; } *ene_vdw = *ene_vdw + evdw; } /*************************************************************************** IMPLEMENTATION OF NB ENERGY Elec Functional NONE Elec CUTOFF NONE VDW Functional LJ VDW CUTOFF SWITCH SOLVATION NONE Sol Functional SolvNone Geometry PLAIN Compute virial no Float type double CPU CPUPlain ***************************************************************************/ void energy_force_hsswn_d (const int *first, const int *second, const int *boundary, const int natoms, const double *coor, double *force, const double c2, const double c_on2, const int *type_code, const int type_size, const double *vdw_data, double *ene_vdw, double *ene_cul) { int i; int j; int f; int fpos; int s; int spos; int b1; int b2; int tc1; int tc2; int data_pos; double fracX; double fracY; double fracZ; double dx, dy, dz; double ix, iy, iz; double jx, jy, jz; double tx, ty, tz; double d2; double d2inv; double d6inv; double d12inv; double sh1; double sh2; double sw; double sw1; double frac; double fvdw; double nb6; double nb12; double enb; double A; double evdw = 0; const double c2diffinv = 1.0 / ((c2 - c_on2) * (c2 - c_on2) * (c2 - c_on2)); const double one = 1.0; const double two = 2.0; const double three = 3.0; const double four = 4.0; const double six = 6.0; const double twelve = 12.0; /* Start outer loop over neighborlists */ for (i = 0; (i < natoms); i++) { /* Initialize outer atom */ f = first[i]; fpos = COOR_DIM * f; ix = coor[fpos]; iy = coor[fpos + 1]; iz = coor[fpos + 2]; /* Load limits for loop over neighbors */ b1 = boundary[2 * i]; b2 = boundary[2 * i + 1]; /* Zero forces */ fracX = 0; fracY = 0; fracZ = 0; /* Load pars */ tc1 = type_code[f]; /* Inner loop */ for (j = b1; (j < b2); j++) { /* Load coor second */ s = second[j]; spos = COOR_DIM * s; jx = coor[spos]; jy = coor[spos + 1]; jz = coor[spos + 2]; /* Compute distances */ dx = ix - jx; dy = iy - jy; dz = iz - jz; d2 = dx * dx + dy * dy + dz * dz; /* Compute force */ if (d2 < c2) { tc2 = type_code[s]; data_pos = 4 * (type_size * tc1 + tc2); A = vdw_data[data_pos]; d2inv = one / d2; d6inv = d2inv * d2inv * d2inv; d12inv = d6inv * d6inv; nb12 = A * d12inv; enb = nb12; fvdw = twelve * nb12; frac = fvdw * d2inv; if (d2 > c_on2) { sh1 = c2 - d2; sh2 = (c2 + two * d2 - three * c_on2); sw = sh1 * sh2; sw1 = four * (sh1 - sh2); frac = frac * sw - (A * d12inv) * sw1; sh1 = c2diffinv * sh1; frac = frac * sh1; sw = sw * sh1; enb = enb * sw; } evdw = evdw + enb; tx = frac * dx; ty = frac * dy; tz = frac * dz; fracX = fracX - tx; fracY = fracY - ty; fracZ = fracZ - tz; force[spos] = force[spos] + tx; force[spos + 1] = force[spos + 1] + ty; force[spos + 2] = force[spos + 2] + tz; } } /* Assign to first */ force[fpos] = force[fpos] + fracX; force[fpos + 1] = force[fpos + 1] + fracY; force[fpos + 2] = force[fpos + 2] + fracZ; } *ene_vdw = *ene_vdw + evdw; } /*************************************************************************** IMPLEMENTATION OF NB ENERGY Elec Functional NONE Elec CUTOFF NONE VDW Functional LJ VDW CUTOFF SWITCH SOLVATION NONE Sol Functional SolvNone Geometry PLAIN Compute virial no Float type float CPU CPUPlain ***************************************************************************/ void energy_force_hsswn_f (const int *first, const int *second, const int *boundary, const int natoms, const float *coor, float *force, const float c2, const float c_on2, const int *type_code, const int type_size, const float *vdw_data, float *ene_vdw, float *ene_cul) { int i; int j; int f; int fpos; int s; int spos; int b1; int b2; int tc1; int tc2; int data_pos; float fracX; float fracY; float fracZ; float dx, dy, dz; float ix, iy, iz; float jx, jy, jz; float tx, ty, tz; float d2; float d2inv; float d6inv; float d12inv; float sh1; float sh2; float sw; float sw1; float frac; float fvdw; float nb6; float nb12; float enb; float A; float evdw = 0; const float c2diffinv = 1.0 / ((c2 - c_on2) * (c2 - c_on2) * (c2 - c_on2)); const float one = 1.0; const float two = 2.0; const float three = 3.0; const float four = 4.0; const float six = 6.0; const float twelve = 12.0; /* Start outer loop over neighborlists */ for (i = 0; (i < natoms); i++) { /* Initialize outer atom */ f = first[i]; fpos = COOR_DIM * f; ix = coor[fpos]; iy = coor[fpos + 1]; iz = coor[fpos + 2]; /* Load limits for loop over neighbors */ b1 = boundary[2 * i]; b2 = boundary[2 * i + 1]; /* Zero forces */ fracX = 0; fracY = 0; fracZ = 0; /* Load pars */ tc1 = type_code[f]; /* Inner loop */ for (j = b1; (j < b2); j++) { /* Load coor second */ s = second[j]; spos = COOR_DIM * s; jx = coor[spos]; jy = coor[spos + 1]; jz = coor[spos + 2]; /* Compute distances */ dx = ix - jx; dy = iy - jy; dz = iz - jz; d2 = dx * dx + dy * dy + dz * dz; /* Compute force */ if (d2 < c2) { tc2 = type_code[s]; data_pos = 4 * (type_size * tc1 + tc2); A = vdw_data[data_pos]; d2inv = one / d2; d6inv = d2inv * d2inv * d2inv; d12inv = d6inv * d6inv; nb12 = A * d12inv; enb = nb12; fvdw = twelve * nb12; frac = fvdw * d2inv; if (d2 > c_on2) { sh1 = c2 - d2; sh2 = (c2 + two * d2 - three * c_on2); sw = sh1 * sh2; sw1 = four * (sh1 - sh2); frac = frac * sw - (A * d12inv) * sw1; sh1 = c2diffinv * sh1; frac = frac * sh1; sw = sw * sh1; enb = enb * sw; } evdw = evdw + enb; tx = frac * dx; ty = frac * dy; tz = frac * dz; fracX = fracX - tx; fracY = fracY - ty; fracZ = fracZ - tz; force[spos] = force[spos] + tx; force[spos + 1] = force[spos + 1] + ty; force[spos + 2] = force[spos + 2] + tz; } } /* Assign to first */ force[fpos] = force[fpos] + fracX; force[fpos + 1] = force[fpos + 1] + fracY; force[fpos + 2] = force[fpos + 2] + fracZ; } *ene_vdw = *ene_vdw + evdw; }