kforce_CDLJ.br

/* -------------------------------------------------------------------------- *
 *                                   OpenMM                                   *
 * -------------------------------------------------------------------------- *
 * This is part of the OpenMM molecular simulation toolkit originating from   *
 * Simbios, the NIH National Center for Physics-Based Simulation of           *
 * Biological Structures at Stanford, funded under the NIH Roadmap for        *
 * Medical Research, grant U54 GM072970. See https://simtk.org.               *
 *                                                                            *
 * Portions copyright (c) 2008 Stanford University and the Authors.           *
 * Authors: Peter Eastman, Mark Friedrichs, Chris Bruns                       *
 * Contributors:                                                              *
 *                                                                            *
 * Permission is hereby granted, free of charge, to any person obtaining a    *
 * copy of this software and associated documentation files (the "Software"), *
 * to deal in the Software without restriction, including without limitation  *
 * the rights to use, copy, modify, merge, publish, distribute, sublicense,   *
 * and/or sell copies of the Software, and to permit persons to whom the      *
 * Software is furnished to do so, subject to the following conditions:       *
 *                                                                            *
 * The above copyright notice and this permission notice shall be included in *
 * all copies or substantial portions of the Software.                        *
 *                                                                            *
 * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR *
 * IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,   *
 * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL    *
 * THE AUTHORS, CONTRIBUTORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM,    *
 * DAMAGES OR OTHER LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR      *
 * OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE  *
 * USE OR OTHER DEALINGS IN THE SOFTWARE.                                     *
 * -------------------------------------------------------------------------- */

kernel float4 scalar_force_CDLJ( float4 qq, float epsfac, float4 sig, float4 eps, float4 r2, float4 params ) {
 float4 invr, invrsig2, invrsig6;
 float4 f;

 invr = rsqrt( r2 );
 invrsig2 = invr * sig;
 invrsig2 = invrsig2 * invrsig2;
 invrsig6 = invrsig2 * invrsig2 * invrsig2;
 f = eps * ( 12.0f*invrsig6 - 6.0f ) * invrsig6;
 f += epsfac*qq*invr;
 f *= invr*invr;
 return f;
}

kernel float scalar_force_single_CDLJ( float qq, float epsfac, float sig, float eps, float r2, float4 params ) {
 float invr, invrsig2, invrsig6;
 float f;

 invr = rsqrt( r2 );
 invrsig2 = invr * sig;
 invrsig2 = invrsig2 * invrsig2;
 invrsig6 = invrsig2 * invrsig2 * invrsig2;
 f = eps * ( 12.0f*invrsig6 - 6.0f ) * invrsig6;
 f += epsfac*qq*invr;
 f *= invr*invr;
 return f;
}
kernel float4 get_r2_CDLJ( float3 d1, float3 d2, float3 d3, float3 d4 ) {
 float4 r2;
 r2 = float4( dot(d1, d1), dot( d2, d2 ), dot( d3, d3 ), dot( d4, d4 ) );
 return r2;
}

kernel void knbforce_CDLJ (
  float natoms,
  float dupfac,
  float strheight,
  float strwidth,
  float jstrwidth,
  float fstrwidth,
  float epsfac,
  float4 params,
  float4 posq[][],
  float4 isigeps<>,
  float4 sigma[][],
  float4 epsilon[][],
  float2 excl[][],


  out float3 outforce1<>,
  out float3 outforce2<> ) {

 float4 jsig, jeps;
 float3 ipos1, ipos2;
 float3 jpos1, jpos2, jpos3, jpos4;
 float qi1, qi2;
 float4 qj, qq;
 float3 d1, d2, d3, d4;
 float2 exclind;
 float4 sig, eps;
 float2 exclusions;
 float4 r2;
 float3 pad;
 float4 exclmask;
 float4 exclconst;
 float2 iatom;
 float2 jstrindex;
 float a_iatom;
 float linind;

 float breakflag;
 float2 jatom;
 float4 fs;
 float jend, jstart;
 float which_rep;

 exclconst = float4( 2.0f, 3.0f, 5.0f, 7.0f );


 outforce1 = float4( 0.0f, 0.0f, 0.0f, 0.0f );
 outforce2 = float4( 0.0f, 0.0f, 0.0f, 0.0f );


   a_iatom = 2.0f*((indexof outforce1).x + (indexof outforce1).y * fstrwidth);


 linind = fmod( a_iatom, natoms );
 iatom.x = fmod( linind, strwidth );
 iatom.y = round( (linind - fmod(linind, strwidth))/strwidth );



 ipos1 = posq[ iatom ].xyz;
 qi1 = posq[ iatom ].w;

 iatom.x += 1;
 ipos2 = posq[ iatom ].xyz;
 qi2 = posq[ iatom ].w;


 breakflag = 1.0f;
 jatom.y = 0.0f;







 which_rep = round( (a_iatom - fmod(a_iatom, natoms))/natoms );
 jend = 1 + round( (natoms - fmod(natoms, (dupfac * strwidth )))/(dupfac * strwidth ) );
 jstart = which_rep * jend;
 if ( which_rep > dupfac - 1 - 0.5f ) {

  jend = 1e6f;
 }
 jend += jstart;

 exclind.x = linind/2.0f;
 exclind.y = jstart * strwidth/4.0f;

 jatom.y = jstart;
 while ( jatom.y < jend && breakflag > 0.0f ) {
  jatom.x = 0;
  while ( jatom.x < strwidth && breakflag > 0.0f ) {

   exclusions = excl[ exclind ];

   exclmask = fmod( exclusions.x, exclconst );


   linind = (jatom.x + jatom.y*strwidth)/4.0f;
   jstrindex.y = round( (linind - fmod(linind, jstrwidth))/jstrwidth );
   jstrindex.x = linind - jstrindex.y * jstrwidth;


   jsig = sigma[ jstrindex ];
   jeps = epsilon[ jstrindex ];

   jpos1 = posq[ jatom ].xyz;
   qj.x = posq[ jatom ].w;
   jatom.x += 1.0f;

   jpos2 = posq[ jatom ].xyz;
   qj.y = posq[ jatom ].w;
   jatom.x += 1.0f;

   jpos3 = posq[ jatom ].xyz;
   qj.z = posq[ jatom ].w;
   jatom.x += 1.0f;

   jpos4 = posq[ jatom ].xyz;
   qj.w = posq[ jatom ].w;
   jatom.x += 1.0f;


   sig = isigeps.x + jsig;
   eps = isigeps.y * jeps;
   qq = qi1 * qj;

   d1 = ipos1 - jpos1;
   d2 = ipos1 - jpos2;
   d3 = ipos1 - jpos3;
   d4 = ipos1 - jpos4;


   r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask * 10000.0f;
   fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );

   outforce1 += fs.x * d1;
   outforce1 += fs.y * d2;
   outforce1 += fs.z * d3;
   outforce1 += fs.w * d4;
/*
outforce1.x += exclmask.x;
outforce1.x += exclmask.y;
outforce1.x += exclmask.z;
outforce1.x += exclmask.w;
if( exclmask.x > 0.5f ){
   outforce1.x += 1.0f;
}
if( exclmask.y > 0.5f ){
   outforce1.x += 1.0f;
}
if( exclmask.z > 0.5f ){
   outforce1.x += 1.0f;
}
if( exclmask.w > 0.5f ){
   outforce1.x += 1.0f;
}
*/

   exclmask = fmod( exclusions.y, exclconst );

   sig = isigeps.z + jsig;
   eps = isigeps.w * jeps;
   qq = qi2 * qj;

   d1 = ipos2 - jpos1;
   d2 = ipos2 - jpos2;
   d3 = ipos2 - jpos3;
   d4 = ipos2 - jpos4;

   r2 = get_r2_CDLJ( d1, d2, d3, d4 ) + exclmask * 10000.0f;
   fs = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );

   outforce2 += fs.x * d1;
   outforce2 += fs.y * d2;
   outforce2 += fs.z * d3;
   outforce2 += fs.w * d4;

/*
outforce2.x += exclmask.x;
outforce2.x += exclmask.y;
outforce2.x += exclmask.z;
outforce2.x += exclmask.w;
if( exclmask.x > 0.5f ){
   outforce2.x += 1.0f;
}
if( exclmask.y > 0.5f ){
   outforce2.x += 1.0f;
}
if( exclmask.z > 0.5f ){
   outforce2.x += 1.0f;
}
if( exclmask.w > 0.5f ){
   outforce2.x += 1.0f;
}
*/

   exclind.y += 1.0f;


   if ( natoms - (jatom.y*strwidth + jatom.x) < 0.9f )
    breakflag = -1.0f;
  }

  jatom.y += 1.0f;
 }

}

kernel void knbforce_CDLJ4(
  float natoms,
  float nAtomsCeiling,
  float dupfac,
  float strheight,
  float strwidth,
  float jstrwidth,
  float fstrwidth,
  float epsfac,
  float4 params,
  float3 posq[][],
  float  charge[][],
  float4 isigeps[][],
  float4 sigma[][],
  float4 epsilon[][],
  float excl[][],
  out float3 outforce1<>,
  out float3 outforce2<>,
  out float3 outforce3<>,
  out float3 outforce4<> ) {

 float4 jsig, jeps;
 float3 ipos1, ipos2, ipos3, ipos4;
 float4 jpos;
 float3 jpos1, jpos2, jpos3, jpos4;
 float qi1, qi2, qi3, qi4;

 float isig1, isig2, isig3, isig4;
 float ieps1, ieps2, ieps3, ieps4;

 float4 qj, qq;
 float3 d1, d2, d3, d4;
 float2 exclind;
 float4 sig, eps;
 float  exclusions;
 float4 r2;
 float3 pad;
 float4 exclmask;
 float4 exclconst;
 float2 iatom;
 float2 jstrindex;
 float a_iatom;
 float linind;

 float breakflag;
 float2 jatom;
 float4 fs;
 float jend, jstart;
 float which_rep;
 float maskMultiplier;

 maskMultiplier = 1.0e+06f;

 exclconst  = float4( 2.0f, 3.0f, 5.0f, 7.0f );

 outforce1  = float3( 0.0f, 0.0f, 0.0f );
 outforce2  = float3( 0.0f, 0.0f, 0.0f );
 outforce3  = float3( 0.0f, 0.0f, 0.0f );
 outforce4  = float3( 0.0f, 0.0f, 0.0f );

 a_iatom    = 4.0f*((indexof outforce1).x + (indexof outforce1).y * fstrwidth);

 // linind     = fmod( a_iatom, natoms );
 linind     = fmod( a_iatom, nAtomsCeiling );

 iatom.x    = fmod( linind, strwidth );
 iatom.y    = round( (linind - fmod(linind, strwidth))/strwidth );

/* --------------------------------------------------------------------- */

 ipos1      = posq[ iatom ];
 qi1        = charge[ iatom ];

 jstrindex  = isigeps[ iatom ];
 isig1      = jstrindex.x;
 ieps1      = jstrindex.y;

/* --------------------------------------------------------------------- */

 iatom.x   += 1;
 ipos2      = posq[ iatom ];
 qi2        = charge[ iatom ];

 jstrindex  = isigeps[ iatom ];
 isig2      = jstrindex.x;
 ieps2      = jstrindex.y;

/* --------------------------------------------------------------------- */

 iatom.x   += 1;
 ipos3      = posq[ iatom ];
 qi3        = charge[ iatom ];

 jstrindex  = isigeps[ iatom ];
 isig3      = jstrindex.x;
 ieps3      = jstrindex.y;

/* --------------------------------------------------------------------- */

 iatom.x   += 1;
 ipos4      = posq[ iatom ];
 qi4        = charge[ iatom ];


 jstrindex  = isigeps[ iatom ];
 isig4      = jstrindex.x;
 ieps4      = jstrindex.y;

/* --------------------------------------------------------------------- */

 breakflag  = 1.0f;
 jatom.y    = 0.0f;

 // which_rep  = round( (a_iatom - fmod(a_iatom, natoms))/natoms );
 which_rep  = round( (a_iatom - fmod(a_iatom, nAtomsCeiling ))/nAtomsCeiling );

 jend       = 1.0f + round( (natoms - fmod(natoms, (dupfac * strwidth )))/(dupfac * strwidth ) );
 jstart     = which_rep * jend;
 if ( which_rep > dupfac - 1.5f ) {
  jend = 1e6f;
 }
 jend      += jstart;

 exclind.x  = linind;
 exclind.y  = jstart * strwidth/4.0f;
 exclind.y  = floor(  exclind.y  + 0.25f );

 jatom.y    = jstart;
 while ( jatom.y < jend && breakflag > 0.0f ) {
  jatom.x = 0;
  while ( jatom.x < strwidth && breakflag > 0.0f ) {

   exclusions  = excl[ exclind ];
   exclmask    = fmod( exclusions, exclconst );

   linind      = (jatom.x + jatom.y*strwidth)/4.0f;
   linind      = floor( linind + 0.25f );

   jstrindex.y = round( (linind - fmod(linind, jstrwidth))/jstrwidth );
   jstrindex.x = linind - jstrindex.y * jstrwidth;
   jstrindex.x = floor( jstrindex.x + 0.25f );

   /* --------------------------------------------------------------------- */

   jsig        = sigma[ jstrindex ];
   jeps        = epsilon[ jstrindex ];

   jpos1       = posq[ jatom ];
   qj.x        = charge[ jatom ];
   jatom.x     = floor( jatom.x + 1.2f );

   jpos2       = posq[ jatom ];
   qj.y        = charge[ jatom ];
   jatom.x     = floor( jatom.x + 1.2f );

   jpos3       = posq[ jatom ];
   qj.z        = charge[ jatom ];
   jatom.x     = floor( jatom.x + 1.2f );

   jpos4       = posq[ jatom ];
   qj.w        = charge[ jatom ];
   jatom.x     = floor( jatom.x + 1.2f );

   /* --------------------------------------------------------------------- */

   sig         = isig1 + jsig;
   eps         = ieps1 * jeps;
   qq          = qi1*qj;

   d1          = ipos1 - jpos1;
   d2          = ipos1 - jpos2;
   d3          = ipos1 - jpos3;
   d4          = ipos1 - jpos4;

   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier);
   fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );

   outforce1  += fs.x * d1;
   outforce1  += fs.y * d2;
   outforce1  += fs.z * d3;
   outforce1  += fs.w * d4;

   /* --------------------------------------------------------------------- */

   //exclind.x  += 1.0f;
   exclind.x   = floor( exclind.x + 1.2f );
   exclusions  = excl[ exclind ];
   exclmask    = fmod( exclusions, exclconst );

   sig         = isig2 + jsig;
   eps         = ieps2 * jeps;
   qq          = qi2 * qj;

   d1          = ipos2 - jpos1;
   d2          = ipos2 - jpos2;
   d3          = ipos2 - jpos3;
   d4          = ipos2 - jpos4;

   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier);
   fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );

   outforce2  += fs.x * d1;
   outforce2  += fs.y * d2;
   outforce2  += fs.z * d3;
   outforce2  += fs.w * d4;

   /* --------------------------------------------------------------------- */

   //exclind.x  += 1.0f;
   exclind.x   = floor( exclind.x + 1.2f );
   exclusions  = excl[ exclind ];
   exclmask    = fmod( exclusions, exclconst );

   sig         = isig3 + jsig;
   eps         = ieps3 * jeps;
   qq          = qi3   * qj;

   d1          = ipos3 - jpos1;
   d2          = ipos3 - jpos2;
   d3          = ipos3 - jpos3;
   d4          = ipos3 - jpos4;

   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier);
   fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );

   outforce3  += fs.x * d1;
   outforce3  += fs.y * d2;
   outforce3  += fs.z * d3;
   outforce3  += fs.w * d4;

   /* --------------------------------------------------------------------- */

   //exclind.x  += 1.0f;
   exclind.x   = floor( exclind.x + 1.2f );
   exclusions  = excl[ exclind ];
   exclmask    = fmod( exclusions, exclconst );

   sig         = isig4 + jsig;
   eps         = ieps4 * jeps;
   qq          = qi4   * qj;

   d1          = ipos4 - jpos1;
   d2          = ipos4 - jpos2;
   d3          = ipos4 - jpos3;
   d4          = ipos4 - jpos4;

   r2          = get_r2_CDLJ( d1, d2, d3, d4 ) + (exclmask*maskMultiplier);
   fs          = scalar_force_CDLJ( qq, epsfac, sig, eps, r2, params );

   outforce4  += fs.x * d1;
   outforce4  += fs.y * d2;
   outforce4  += fs.z * d3;
   outforce4  += fs.w * d4;

   /* --------------------------------------------------------------------- */

   //exclind.x  -= 3.0f;
   exclind.x   = floor( exclind.x - 2.9f );

   //exclind.y  += 1.0f;
   exclind.y  = floor(  exclind.y  + 1.25f );

   if( natoms - (jatom.y*strwidth + jatom.x) < 0.9f )breakflag = -1.0f;
  }

  //jatom.y += 1.0f;
  jatom.y  = floor( jatom.y + 1.2f );
 }

}

kernel void kforce14_CDLJ(
  float natoms,
  float strwidth,
  float epsfac,
  float4 params,
  float2 atoms<>,
  float4 posq[][],
  float2 sigeps<>,
  out float3 fi<>,
  out float3 fj<>
  )
{
 float2 ai, aj;
 float3 d1;
 float r2;
 float qq, fs;

 //ai.y = floor( atoms.x / strwidth );
 ai.y = round( (atoms.x - fmod( atoms.x, strwidth ))/strwidth );
 ai.x = atoms.x - ai.y * strwidth;

 //aj.y = floor( atoms.y / strwidth );
 aj.y = round( (atoms.y - fmod( atoms.y, strwidth ))/strwidth );
 aj.x = atoms.y - aj.y * strwidth;

 d1 = posq[ai].xyz - posq[aj].xyz;
 qq = posq[ai].w * posq[aj].w;

 r2 = dot( d1, d1 );

 fs = scalar_force_single_CDLJ( qq, epsfac, sigeps.x, sigeps.y, r2, params );

 fi = fs.x * d1;
 fj = -fi;
}

kernel void kbonded_CDLJ (
  float epsfac,
  float xstrwidth,
  float4 nbparams,
  float3 posq[][],
  float  charge[][],
  float4 atoms<>,
  float4 parm0<>,
  float4 parm1<>,
  float4 parm2<>,
  float4 parm3<>,
  float4 parm4<>,
  out float3 fi<>,
  out float3 fj<>,
  out float3 fk<>,
  out float3 fl<>
  )
{
 float3 rij, rkj, rkl, ril;
 float2 ai, aj, ak, al;
 float3 m, n;
 float sgnphi;
 float cosfac;
 float ddphi, ddphi_rb;
 float3 u, v, s;
 float invlkj, invlij, invlkl, msq, nsq, cos_phi, sin_phi;
 float3 uij, ukj, ukl;
 float phi, mdphi;
 float rij_d_ukj, ukj_d_rkl;
 float normfac;
 float qq;
 float3 fi_ang, fj_ang, fk_ang;
 float3 fi_bond;
 float cii, ckk, cik;
 float fs, st, sth;
 float3 fi_pair;
 float r2 ;
 float torsion_numerator;
 float theta1, costheta1, invsintheta1;
 float theta2, costheta2, invsintheta2;
 float3 posqi, posqj, posqk, posql;
 float  chargei, chargel;
 float sig, eps;






   ai.y = round( (atoms.x - fmod( atoms.x, xstrwidth ))/xstrwidth );
 ai.x = atoms.x - ai.y * xstrwidth;
 posqi = posq[ ai ];
 chargei = charge[ ai ];






 if ( atoms.y > -0.5f ) {

      aj.y = round( (atoms.y - fmod( atoms.y, xstrwidth ))/xstrwidth );
  aj.x = atoms.y - aj.y * xstrwidth;
  posqj = posq[ aj ];
 }
 else
  posqj = float4( 0.0f, 0.0f, 1.0f, 0.0f );

 if ( atoms.z > -0.5f ) {

      ak.y = round( (atoms.z - fmod( atoms.z, xstrwidth ))/xstrwidth );
  ak.x = atoms.z - ak.y * xstrwidth;
  posqk = posq[ ak ];
 }
 else
  posqk = float4( 0.0f, 1.0f, 1.0f, 0.0f );

 if ( atoms.w > -0.5f ) {

      al.y = round( (atoms.w - fmod( atoms.w, xstrwidth ))/xstrwidth );
  al.x = atoms.w - al.y * xstrwidth;
  posql = posq[ al ];
  chargel = charge[ al ];
 } else {
  posql = float4( 1.0f, 1.0f, 1.0f, 0.0f );
  chargel = 0.0f;
}

 qq = chargei*chargel;
 rij = posqi.xyz - posqj.xyz;
 rkj = posqk.xyz - posqj.xyz;
 rkl = posqk.xyz - posql.xyz;
 ril = posqi.xyz - posql.xyz;



 invlij = rsqrt( dot( rij, rij ) );
 invlkj = rsqrt( dot( rkj, rkj ) );
 invlkl = rsqrt( dot( rkl, rkl ) );

 uij = rij * invlij;
 ukj = rkj * invlkj;
 ukl = rkl * invlkl;


 rij_d_ukj = dot( rij, ukj );
 ukj_d_rkl = dot( ukj, rkl );

 m = cross( uij, ukj );
 n = cross( ukj, ukl );


 costheta1 = clamp( rij_d_ukj * invlij, -1.0f, 1.0f );
 theta1 = acos( costheta1 );
 invsintheta1 = rsqrt( 1.0f - costheta1 * costheta1 );


 costheta2 = clamp( ukj_d_rkl * invlkl, -1.0f, 1.0f );
 theta2 = acos( costheta2 );
 invsintheta2 = rsqrt( 1.0f - costheta2 * costheta2 );

 cos_phi = clamp( dot(m,n) * invsintheta1 * invsintheta2, -1.0f, 1.0f );

 sgnphi = sign( dot(rij, n) );


 phi   = sgnphi * acos( cos_phi );
 mdphi = parm1.w * phi - parm1.z;
 ddphi = -parm1.y * parm1.w * sin( mdphi );


 cos_phi = -cos_phi;
 sin_phi = -sgnphi*sqrt( clamp( 1.0f - cos_phi * cos_phi, 0.0f, 1.0f) );

 ddphi_rb = 5.0f * parm1.x;
 ddphi_rb = 4.0f * parm0.w + ddphi_rb * cos_phi;
 ddphi_rb = 3.0f * parm0.z + ddphi_rb * cos_phi;
 ddphi_rb = 2.0f * parm0.y + ddphi_rb * cos_phi;
 ddphi_rb = parm0.x + ddphi_rb * cos_phi;
 ddphi_rb = -ddphi_rb * sin_phi;


 ddphi += ddphi_rb;

 fi = (-ddphi * invlij * invsintheta1 * invsintheta1 ) * m;
 fl = ( ddphi * invlkl * invsintheta2 * invsintheta2 ) * n;


 u = rij_d_ukj * invlkj * fi;
 v = ukj_d_rkl * invlkj * fl;
 s = u - v;

 fj = s - fi;
 fk = -(s + fl);





 fs = -parm2.y * ( theta1 - parm2.x );
 st = fs * invsintheta1;
 st = clamp( st, -100000.0f, 100000.0f );

 sth = st * costheta1;

 cik = st * invlkj * invlij;
 cii = sth * invlij;
 ckk = sth * invlkj;

 fi_ang = -( cik * rkj - cii * uij );
 fk_ang = -( cik * rij - ckk * ukj );
 fj_ang = -fi_ang - fk_ang;

 fi += fi_ang;
 fj += fj_ang;
 fk += fk_ang;
 fs = -parm2.w * ( theta2 - parm2.z );
 st = fs * invsintheta2;
 st = clamp( st, -100000.0f, 100000.0f );

 sth = st * costheta2;

 cik = st * invlkj * invlkl;
 cii = sth * invlkj;
 ckk = sth * invlkl;

 fi_ang = ( cik * rkl - cii * ukj );
 fk_ang = ( cik * rkj - ckk * ukl );
 fj_ang = -fi_ang - fk_ang;

 fj += fi_ang;
 fk += fj_ang;
 fl += fk_ang;


 fi_bond = -parm3.y * ( 1.0f - parm3.x * invlij ) * rij;
 fi += fi_bond;
 fj += -fi_bond;


 fi_bond = parm3.w * ( 1.0f - parm3.z * invlkj ) * rkj;
 fj += fi_bond;
 fk += -fi_bond;


 fi_bond = -parm4.y * ( 1.0f - parm4.x * invlkl ) * rkl;
 fk += fi_bond;
 fl += -fi_bond;

 if ( parm4.z > -0.5f ) {
  r2 = dot( ril, ril );
  fs = scalar_force_single_CDLJ( qq, epsfac, parm4.z, parm4.w, r2, nbparams );
  fi_pair = fs * ril;
  fi += fi_pair;
  fl -= fi_pair;
 }

}