/*************************************************************************** forcefuncsoft.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_softswn_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 double alpha, 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 sh1; double sh2; double sw; double sw1; double frac; double nbvdw; double A; double B; int tc1,tc2; int data_pos; 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_inv = 1.0/6.0; const double twelve_inv = 1.0/12.0; int i,j; // printf("DEBUG SOFT %f\n",alpha); 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-six_inv*B*d6inv)*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 energy_force_softswn_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 double alpha, 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 B; 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; // printf("DEBUG SOFT %f\n",alpha); /* 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]; B = vdw_data[data_pos + 1]; d2inv = one / (d2+alpha); d6inv = d2inv * d2inv * d2inv; d12inv = d6inv * d6inv; nb12 = A * d12inv; nb6 = B * d6inv; enb = nb12-nb6; fvdw = twelve * nb12-six*nb6; 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- B * d6inv) * 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; } //float void force_softswn_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 float alpha, 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; float B; 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-six_inv*B*d6inv)*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 energy_force_softswn_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 float alpha, 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 B; 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]; B = vdw_data[data_pos + 1]; d2inv = one / (d2+alpha); d6inv = d2inv * d2inv * d2inv; d12inv = d6inv * d6inv; nb12 = A * d12inv; nb6 = B * d6inv; enb = nb12-nb6; fvdw = twelve * nb12-six*nb6; 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- B * d6inv) * 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; }