/***************************************************************************
                          pbcenefunc.cpp  -  description
                             -------------------
    begin                : Mon Nov 15 23:15:24 2004
    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

#include <math.h>
#include <almost.h>
#ifdef IX86
#define x86trunc(a,b) asm("fld %1 ; fistpl %0" : "=m" (*&b) : "f" (a));
#else
#define x86trunc(a,b) {b = (int)a;}
#endif

//cssn
void pbcenergy_cssn_d(const int * first,
		      const int * second,
		      const int * boundarry,
		      const int natoms,
		      const double * coor,
		      const double  c2,
		      const double * charge,
		      const int * type_code,
		      const int type_size,
		      const double * vdw_data,
		      double * ene_vdw,
		      double * ene_cul){

  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double jx,jy,jz;
  double d2;
  double d6;
  double dinv;
  double d2inv;
  double d6inv;
  double d12inv;
  double shiftf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double c2inv = 1.0/c2;
  const double c6inv = c2inv*c2inv*c2inv;
  const double c12inv = c6inv*c6inv;
  const double one = 1.0;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      //	  if(d2>c2){ cout<<"AHHHHHHHHHH"<<endl;}
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	dinv = sqrt(d2inv);
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*dinv;
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void pbcenergy_cssn_f(const int * first,
		      const int * second,
		      const int * boundarry,
		      const int natoms,
		      const float * coor,
		      const float  c2,
		      const float * charge,
		      const int * type_code,
		      const int type_size,
		      const float * vdw_data,
		      float * ene_vdw,
		      float * ene_cul){

  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float jx,jy,jz;
  float d2;
  float d6;
  float dinv;
  float d2inv;
  float d6inv;
  float d12inv;
  float shiftf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float c2inv = 1.0/c2;
  const float c6inv = c2inv*c2inv*c2inv;
  const float c12inv = c6inv*c6inv;
  const float one = 1.0;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      //	  if(d2>c2){ cout<<"AHHHHHHHHHH"<<endl;}
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	dinv = sqrt(d2inv);
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*dinv;
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

void pbcenergy_sh_cssn_d(const int * first,
			 const int * second,
			 const int * boundarry,
			 const int natoms,
			 const double * coor,
			 const double * sh,
			 const double  c2,
			 const double * charge,
			 const int * type_code,
			 const int type_size,
			 const double * vdw_data,
			 double * ene_vdw,
			 double * ene_cul){

  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double shX,shY,shZ;
  double jx,jy,jz;
  double d2;
  double d6;
  double dinv;
  double d2inv;
  double d6inv;
  double d12inv;
  double shiftf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double c2inv = 1.0/c2;
  const double c6inv = c2inv*c2inv*c2inv;
  const double c12inv = c6inv*c6inv;
  const double one = 1.0;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      //	  if(d2>c2){ cout<<"AHHHHHHHHHH"<<endl;}
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	dinv = sqrt(d2inv);
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*dinv;
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void pbcenergy_sh_cssn_f(const int * first,
			 const int * second,
			 const int * boundarry,
			 const int natoms,
			 const float * coor,
			 const float * sh,
			 const float  c2,
			 const float * charge,
			 const int * type_code,
			 const int type_size,
			 const float * vdw_data,
			 float * ene_vdw,
			 float * ene_cul){

  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float shX,shY,shZ;
  float jx,jy,jz;
  float d2;
  float d6;
  float dinv;
  float d2inv;
  float d6inv;
  float d12inv;
  float shiftf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float c2inv = 1.0/c2;
  const float c6inv = c2inv*c2inv*c2inv;
  const float c12inv = c6inv*c6inv;
  const float one = 1.0;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      //	  if(d2>c2){ cout<<"AHHHHHHHHHH"<<endl;}
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	dinv = sqrt(d2inv);
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*dinv;
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}



//cswswn
void 
pbcenergy_cswswn_d(const int * first,
		   const int * second,
		   const int * boundarry,
		   const int natoms,
		   const double * coor,
		   const double  c2,
		   const double  c_on2,
		   const double * charge,
		   const int * type_code,
		   const int type_size,
		   const double * vdw_data,
		   double * ene_vdw,
		   double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double jx,jy,jz;
  double d2;
  double dinv;
  double d2inv;
  double d6inv;
  double d12inv;
  double switchf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	dinv = sqrt(d2inv);
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;
	  ecul = ecul + q1*q2*switchf*dinv;
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	  ecul = ecul + q1*q2*dinv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void 
pbcenergy_cswswn_f(const int * first,
		   const int * second,
		   const int * boundarry,
		   const int natoms,
		   const float * coor,
		   const float  c2,
		   const float  c_on2,
		   const float * charge,
		   const int * type_code,
		   const int type_size,
		   const float * vdw_data,
		   float * ene_vdw,
		   float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float jx,jy,jz;
  float d2;
  float dinv;
  float d2inv;
  float d6inv;
  float d12inv;
  float switchf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	dinv = sqrt(d2inv);
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;
	  ecul = ecul + q1*q2*switchf*dinv;
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	  ecul = ecul + q1*q2*dinv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void 
pbcenergy_sh_cswswn_d(const int * first,
		      const int * second,
		      const int * boundarry,
		      const int natoms,
		      const double * coor,
		      const double * sh,
		      const double  c2,
		      const double  c_on2,
		      const double * charge,
		      const int * type_code,
		      const int type_size,
		      const double * vdw_data,
		      double * ene_vdw,
		      double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double shX,shY,shZ;
  double jx,jy,jz;
  double d2;
  double dinv;
  double d2inv;
  double d6inv;
  double d12inv;
  double switchf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	dinv = sqrt(d2inv);
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;
	  ecul = ecul + q1*q2*switchf*dinv;
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	  ecul = ecul + q1*q2*dinv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void 
pbcenergy_sh_cswswn_f(const int * first,
		      const int * second,
		      const int * boundarry,
		      const int natoms,
		      const float * coor,
		      const float * sh,
		      const float  c2,
		      const float  c_on2,
		      const float * charge,
		      const int * type_code,
		      const int type_size,
		      const float * vdw_data,
		      float * ene_vdw,
		      float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float shX,shY,shZ;
  float jx,jy,jz;
  float d2;
  float dinv;
  float d2inv;
  float d6inv;
  float d12inv;
  float switchf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	dinv = sqrt(d2inv);
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;
	  ecul = ecul + q1*q2*switchf*dinv;
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	  ecul = ecul + q1*q2*dinv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

//csswn
void 
pbcenergy_csswn_d(const int * first,
		  const int * second,
		  const int * boundarry,
		  const int natoms,
		  const double * coor,
		  const double  c2,
		  const double  c_on2,
		  const double * charge,
		  const int * type_code,
		  const int type_size,
		  const double * vdw_data,
		  double * ene_vdw,
		  double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double jx,jy,jz;
  double d2;
  double dinv;
  double d2inv;
  double d6inv;
  double d12inv;
  double shiftf;
  double switchf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double c2inv = 1.0/c2;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	dinv = sqrt(d2inv);
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*dinv;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;	  
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void 
pbcenergy_csswn_f(const int * first,
		  const int * second,
		  const int * boundarry,
		  const int natoms,
		  const float * coor,
		  const float  c2,
		  const float  c_on2,
		  const float * charge,
		  const int * type_code,
		  const int type_size,
		  const float * vdw_data,
		  float * ene_vdw,
		  float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float jx,jy,jz;
  float d2;
  float dinv;
  float d2inv;
  float d6inv;
  float d12inv;
  float shiftf;
  float switchf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float c2inv = 1.0/c2;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	dinv = sqrt(d2inv);
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*dinv;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;	  
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

void 
pbcenergy_sh_csswn_d(const int * first,
		     const int * second,
		     const int * boundarry,
		     const int natoms,
		     const double * coor,
		     const double * sh,
		     const double  c2,
		     const double  c_on2,
		     const double * charge,
		     const int * type_code,
		     const int type_size,
		     const double * vdw_data,
		     double * ene_vdw,
		     double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double shX,shY,shZ;
  double jx,jy,jz;
  double d2;
  double dinv;
  double d2inv;
  double d6inv;
  double d12inv;
  double shiftf;
  double switchf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double c2inv = 1.0/c2;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	dinv = sqrt(d2inv);
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*dinv;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;	  
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void 
pbcenergy_sh_csswn_f(const int * first,
		     const int * second,
		     const int * boundarry,
		     const int natoms,
		     const float * coor,
		     const float * sh,
		     const float  c2,
		     const float  c_on2,
		     const float * charge,
		     const int * type_code,
		     const int type_size,
		     const float * vdw_data,
		     float * ene_vdw,
		     float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float shX,shY,shZ;
  float jx,jy,jz;
  float d2;
  float dinv;
  float d2inv;
  float d6inv;
  float d12inv;
  float shiftf;
  float switchf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float c2inv = 1.0/c2;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	dinv = sqrt(d2inv);
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*dinv;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;	  
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

//cswsn
void 
pbcenergy_cswsn_d(const int * first,
		  const int * second,
		  const int * boundarry,
		  const int natoms,
		  const double * coor,
		  const double  c2,
		  const double  c_on2,
		  const double * charge,
		  const int * type_code,
		  const int type_size,
		  const double * vdw_data,
		  double * ene_vdw,
		  double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double jx,jy,jz;
  double d2;
  double d6;
  double dinv;
  double d2inv;
  double d6inv;
  double d12inv;
  double shiftf;
  double switchf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double c2inv = 1.0/c2;
  const double c6inv = c2inv*c2inv*c2inv;
  const double c12inv = c6inv*c6inv;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	dinv = sqrt(d2inv);
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  ecul = ecul + q1*q2*switchf*dinv;
	}
	else {
	  ecul = ecul + q1*q2*dinv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

void 
pbcenergy_cswsn_f(const int * first,
		  const int * second,
		  const int * boundarry,
		  const int natoms,
		  const float * coor,
		  const float  c2,
		  const float  c_on2,
		  const float * charge,
		  const int * type_code,
		  const int type_size,
		  const float * vdw_data,
		  float * ene_vdw,
		  float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float jx,jy,jz;
  float d2;
  float d6;
  float dinv;
  float d2inv;
  float d6inv;
  float d12inv;
  float shiftf;
  float switchf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float c2inv = 1.0/c2;
  const float c6inv = c2inv*c2inv*c2inv;
  const float c12inv = c6inv*c6inv;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	dinv = sqrt(d2inv);
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  ecul = ecul + q1*q2*switchf*dinv;
	}
	else {
	  ecul = ecul + q1*q2*dinv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

void 
pbcenergy_sh_cswsn_d(const int * first,
		     const int * second,
		     const int * boundarry,
		     const int natoms,
		     const double * coor,
		     const double * sh,
		     const double  c2,
		     const double  c_on2,
		     const double * charge,
		     const int * type_code,
		     const int type_size,
		     const double * vdw_data,
		     double * ene_vdw,
		     double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double shX,shY,shZ;
  double jx,jy,jz;
  double d2;
  double d6;
  double dinv;
  double d2inv;
  double d6inv;
  double d12inv;
  double shiftf;
  double switchf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double c2inv = 1.0/c2;
  const double c6inv = c2inv*c2inv*c2inv;
  const double c12inv = c6inv*c6inv;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	dinv = sqrt(d2inv);
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  ecul = ecul + q1*q2*switchf*dinv;
	}
	else {
	  ecul = ecul + q1*q2*dinv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

void 
pbcenergy_sh_cswsn_f(const int * first,
		     const int * second,
		     const int * boundarry,
		     const int natoms,
		     const float * coor,
		     const float * sh,
		     const float  c2,
		     const float  c_on2,
		     const float * charge,
		     const int * type_code,
		     const int type_size,
		     const float * vdw_data,
		     float * ene_vdw,
		     float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float shX,shY,shZ;
  float jx,jy,jz;
  float d2;
  float d6;
  float dinv;
  float d2inv;
  float d6inv;
  float d12inv;
  float shiftf;
  float switchf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float c2inv = 1.0/c2;
  const float c6inv = c2inv*c2inv*c2inv;
  const float c12inv = c6inv*c6inv;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	dinv = sqrt(d2inv);
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  ecul = ecul + q1*q2*switchf*dinv;
	}
	else {
	  ecul = ecul + q1*q2*dinv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


//cttn
void 
pbcenergy_cttn_d(const int * first,
		 const int * second,
		 const int * boundarry,
		 const int natoms,
		 const double * coor,
		 const double  c2,
		 const double * charge,
		 const int * type_code,
		 const int type_size,
		 const double * vdw_data,
		 const double * sdata,
		 double * ene_vdw,
		 double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double jx,jy,jz;
  double d2;
  double d;
  double dinv;
  double dt;
  double A,B,q1,q2;
  double eps,eps2,eps3;
  int tc1,tc2;
  int data_pos;
  int spline_pos;
  int n0;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double one = 1.0;
  const double two = 2.0;
  const double three = 3.0;
  const double density = 2000;
  /*   const double h = (double)1/(double)density; */
  /*   const double hinv = 1/h; */
  /*   const int st = (int)hinv; */
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	dinv = invsqrt(d2);	
	d = d2*dinv;
	dt = density*d;
	x86trunc(dt,n0);
	spline_pos = 12*n0;
	eps = dt-n0;
	eps2 = eps*eps;
	eps3 = eps2*eps;
	ecul = ecul + q1*q2*(sdata[spline_pos] + 
			     eps*sdata[spline_pos+1]+ 
			     eps2*sdata[spline_pos+2] +
			     eps3*sdata[spline_pos+3]); 
	nb6 = B*(sdata[spline_pos+4] + 
		 eps*sdata[spline_pos+5]+ 
		 eps2*sdata[spline_pos+6] +
		 eps3*sdata[spline_pos+7]);
	nb12 = A*(sdata[spline_pos+8] + 
		  eps*sdata[spline_pos+9]+ 
		  eps2*sdata[spline_pos+10] +
		  eps3*sdata[spline_pos+11]);
	evdw = evdw + nb12-nb6;
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void 
pbcenergy_cttn_f(const int * first,
		 const int * second,
		 const int * boundarry,
		 const int natoms,
		 const float * coor,
		 const float c2,
		 const float * charge,
		 const int * type_code,
		 const int type_size,
		 const float * vdw_data,
		 const float * sdata,
		 float * ene_vdw,
		 float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float jx,jy,jz;
  float d2;
  float d;
  float dinv;
  float dt;
  float A,B,q1,q2;
  float eps,eps2,eps3;
  int tc1,tc2;
  int data_pos;
  int spline_pos;
  int n0;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float one = 1.0;
  const float two = 2.0;
  const float three = 3.0;
  const float density = 500;
  /*   const float h = (float)1/(float)density; */
  /*   const float hinv = 1/h; */
  /*   const int st = (int)hinv; */
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	dinv = invsqrt(d2);	
	d = d2*dinv;
	dt = density*d;
	x86trunc(dt,n0);
	//n0 = (int)dt;
	spline_pos = 12*n0;
	eps = dt-n0;
	eps2 = eps*eps;
	eps3 = eps2*eps;
	ecul = ecul + q1*q2*(sdata[spline_pos] + 
			     eps*sdata[spline_pos+1]+ 
			     eps2*sdata[spline_pos+2] +
			     eps3*sdata[spline_pos+3]); 
	spline_pos = spline_pos + 4;
	nb6 = B*(sdata[spline_pos] + 
		 eps*sdata[spline_pos+1]+ 
		 eps2*sdata[spline_pos+2] +
		 eps3*sdata[spline_pos+3]);
	spline_pos = spline_pos + 4;
	nb12 = A*(sdata[spline_pos] + 
		  eps*sdata[spline_pos+1]+ 
		  eps2*sdata[spline_pos+2] +
		  eps3*sdata[spline_pos+3]);
	evdw = evdw + nb12-nb6;
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

void 
pbcenergy_sh_cttn_d(const int * first,
		    const int * second,
		    const int * boundarry,
		    const int natoms,
		    const double * coor,
		    const double * sh,
		    const double  c2,
		    const double * charge,
		    const int * type_code,
		    const int type_size,
		    const double * vdw_data,
		    const double * sdata,
		    double * ene_vdw,
		    double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double shX,shY,shZ;
  double jx,jy,jz;
  double d2;
  double d;
  double dinv;
  double dt;
  double A,B,q1,q2;
  double eps,eps2,eps3;
  int tc1,tc2;
  int data_pos;
  int spline_pos;
  int n0;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double one = 1.0;
  const double two = 2.0;
  const double three = 3.0;
  const double density = 2000;
  /*   const double h = (double)1/(double)density; */
  /*   const double hinv = 1/h; */
  /*   const int st = (int)hinv; */
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	dinv = invsqrt(d2);	
	d = d2*dinv;
	dt = density*d;
	x86trunc(dt,n0);
	spline_pos = 12*n0;
	eps = dt-n0;
	eps2 = eps*eps;
	eps3 = eps2*eps;
	ecul = ecul + q1*q2*(sdata[spline_pos] + 
			     eps*sdata[spline_pos+1]+ 
			     eps2*sdata[spline_pos+2] +
			     eps3*sdata[spline_pos+3]); 
	nb6 = B*(sdata[spline_pos+4] + 
		 eps*sdata[spline_pos+5]+ 
		 eps2*sdata[spline_pos+6] +
		 eps3*sdata[spline_pos+7]);
	nb12 = A*(sdata[spline_pos+8] + 
		  eps*sdata[spline_pos+9]+ 
		  eps2*sdata[spline_pos+10] +
		  eps3*sdata[spline_pos+11]);
	evdw = evdw + nb12-nb6;
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void 
pbcenergy_sh_cttn_f(const int * first,
		    const int * second,
		    const int * boundarry,
		    const int natoms,
		    const float * coor,
		    const float * sh,
		    const float c2,
		    const float * charge,
		    const int * type_code,
		    const int type_size,
		    const float * vdw_data,
		    const float * sdata,
		    float * ene_vdw,
		    float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float shX,shY,shZ;
  float jx,jy,jz;
  float d2;
  float d;
  float dinv;
  float dt;
  float A,B,q1,q2;
  float eps,eps2,eps3;
  int tc1,tc2;
  int data_pos;
  int spline_pos;
  int n0;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float one = 1.0;
  const float two = 2.0;
  const float three = 3.0;
  const float density = 500;
  /*   const float h = (float)1/(float)density; */
  /*   const float hinv = 1/h; */
  /*   const int st = (int)hinv; */
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	dinv = invsqrt(d2);	
	d = d2*dinv;
	dt = density*d;
	x86trunc(dt,n0);
	//n0 = (int)dt;
	spline_pos = 12*n0;
	eps = dt-n0;
	eps2 = eps*eps;
	eps3 = eps2*eps;
	ecul = ecul + q1*q2*(sdata[spline_pos] + 
			     eps*sdata[spline_pos+1]+ 
			     eps2*sdata[spline_pos+2] +
			     eps3*sdata[spline_pos+3]); 
	spline_pos = spline_pos + 4;
	nb6 = B*(sdata[spline_pos] + 
		 eps*sdata[spline_pos+1]+ 
		 eps2*sdata[spline_pos+2] +
		 eps3*sdata[spline_pos+3]);
	spline_pos = spline_pos + 4;
	nb12 = A*(sdata[spline_pos] + 
		  eps*sdata[spline_pos+1]+ 
		  eps2*sdata[spline_pos+2] +
		  eps3*sdata[spline_pos+3]);
	evdw = evdw + nb12-nb6;
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}






//rssn
void pbcenergy_rssn_d(const int * first,
		      const int * second,
		      const int * boundarry,
		      const int natoms,
		      const double * coor,
		      const double  c2,
		      const double * charge,
		      const int * type_code,
		      const int type_size,
		      const double * vdw_data,
		      double * ene_vdw,
		      double * ene_cul){

  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double jx,jy,jz;
  double d2;
  double d6;
  double d2inv;
  double d6inv;
  double d12inv;
  double shiftf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double c2inv = 1.0/c2;
  const double c6inv = c2inv*c2inv*c2inv;
  const double c12inv = c6inv*c6inv;
  const double one = 1.0;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      //	  if(d2>c2){ cout<<"AHHHHHHHHHH"<<endl;}
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*d2inv;
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void pbcenergy_rssn_f(const int * first,
		      const int * second,
		      const int * boundarry,
		      const int natoms,
		      const float * coor,
		      const float  c2,
		      const float * charge,
		      const int * type_code,
		      const int type_size,
		      const float * vdw_data,
		      float * ene_vdw,
		      float * ene_cul){

  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float jx,jy,jz;
  float d2;
  float d6;
  float d2inv;
  float d6inv;
  float d12inv;
  float shiftf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float c2inv = 1.0/c2;
  const float c6inv = c2inv*c2inv*c2inv;
  const float c12inv = c6inv*c6inv;
  const float one = 1.0;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      //	  if(d2>c2){ cout<<"AHHHHHHHHHH"<<endl;}
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*d2inv;
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

void pbcenergy_sh_rssn_d(const int * first,
			 const int * second,
			 const int * boundarry,
			 const int natoms,
			 const double * coor,
			 const double * sh,
			 const double  c2,
			 const double * charge,
			 const int * type_code,
			 const int type_size,
			 const double * vdw_data,
			 double * ene_vdw,
			 double * ene_cul){

  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double shX,shY,shZ;
  double jx,jy,jz;
  double d2;
  double d6;
  double d2inv;
  double d6inv;
  double d12inv;
  double shiftf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double c2inv = 1.0/c2;
  const double c6inv = c2inv*c2inv*c2inv;
  const double c12inv = c6inv*c6inv;
  const double one = 1.0;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      //	  if(d2>c2){ cout<<"AHHHHHHHHHH"<<endl;}
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*d2inv;
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void pbcenergy_sh_rssn_f(const int * first,
			 const int * second,
			 const int * boundarry,
			 const int natoms,
			 const float * coor,
			 const float * sh,
			 const float  c2,
			 const float * charge,
			 const int * type_code,
			 const int type_size,
			 const float * vdw_data,
			 float * ene_vdw,
			 float * ene_cul){

  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float shX,shY,shZ;
  float jx,jy,jz;
  float d2;
  float d6;
  float d2inv;
  float d6inv;
  float d12inv;
  float shiftf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float c2inv = 1.0/c2;
  const float c6inv = c2inv*c2inv*c2inv;
  const float c12inv = c6inv*c6inv;
  const float one = 1.0;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      //	  if(d2>c2){ cout<<"AHHHHHHHHHH"<<endl;}
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*d2inv;
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

//rswswn
void 
pbcenergy_rswswn_d(const int * first,
		   const int * second,
		   const int * boundarry,
		   const int natoms,
		   const double * coor,
		   const double  c2,
		   const double  c_on2,
		   const double * charge,
		   const int * type_code,
		   const int type_size,
		   const double * vdw_data,
		   double * ene_vdw,
		   double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double jx,jy,jz;
  double d2;
  double d2inv;
  double d6inv;
  double d12inv;
  double switchf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;
	  ecul = ecul + q1*q2*switchf*d2inv;
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	  ecul = ecul + q1*q2*d2inv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void 
pbcenergy_rswswn_f(const int * first,
		   const int * second,
		   const int * boundarry,
		   const int natoms,
		   const float * coor,
		   const float  c2,
		   const float  c_on2,
		   const float * charge,
		   const int * type_code,
		   const int type_size,
		   const float * vdw_data,
		   float * ene_vdw,
		   float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float jx,jy,jz;
  float d2;
  float d2inv;
  float d6inv;
  float d12inv;
  float switchf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;
	  ecul = ecul + q1*q2*switchf*d2inv;
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	  ecul = ecul + q1*q2*d2inv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

void 
pbcenergy_sh_rswswn_d(const int * first,
		      const int * second,
		      const int * boundarry,
		      const int natoms,
		      const double * coor,
		      const double * sh,
		      const double  c2,
		      const double  c_on2,
		      const double * charge,
		      const int * type_code,
		      const int type_size,
		      const double * vdw_data,
		      double * ene_vdw,
		      double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double shX,shY,shZ;
  double jx,jy,jz;
  double d2;
  double d2inv;
  double d6inv;
  double d12inv;
  double switchf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;
	  ecul = ecul + q1*q2*switchf*d2inv;
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	  ecul = ecul + q1*q2*d2inv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void 
pbcenergy_sh_rswswn_f(const int * first,
		      const int * second,
		      const int * boundarry,
		      const int natoms,
		      const float * coor,
		      const float * sh,
		      const float  c2,
		      const float  c_on2,
		      const float * charge,
		      const int * type_code,
		      const int type_size,
		      const float * vdw_data,
		      float * ene_vdw,
		      float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float shX,shY,shZ;
  float jx,jy,jz;
  float d2;
  float d2inv;
  float d6inv;
  float d12inv;
  float switchf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;
	  ecul = ecul + q1*q2*switchf*d2inv;
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	  ecul = ecul + q1*q2*d2inv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


//rsswn
void 
pbcenergy_rsswn_d(const int * first,
		  const int * second,
		  const int * boundarry,
		  const int natoms,
		  const double * coor,
		  const double  c2,
		  const double  c_on2,
		  const double * charge,
		  const int * type_code,
		  const int type_size,
		  const double * vdw_data,
		  double * ene_vdw,
		  double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double jx,jy,jz;
  double d2;
  double d2inv;
  double d6inv;
  double d12inv;
  double shiftf;
  double switchf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double c2inv = 1.0/c2;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*d2inv;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;	  
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void 
pbcenergy_rsswn_f(const int * first,
		  const int * second,
		  const int * boundarry,
		  const int natoms,
		  const float * coor,
		  const float  c2,
		  const float  c_on2,
		  const float * charge,
		  const int * type_code,
		  const int type_size,
		  const float * vdw_data,
		  float * ene_vdw,
		  float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float jx,jy,jz;
  float d2;
  float d2inv;
  float d6inv;
  float d12inv;
  float shiftf;
  float switchf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float c2inv = 1.0/c2;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*d2inv;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;	  
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

void 
pbcenergy_sh_rsswn_d(const int * first,
		     const int * second,
		     const int * boundarry,
		     const int natoms,
		     const double * coor,
		     const double * sh,
		     const double  c2,
		     const double  c_on2,
		     const double * charge,
		     const int * type_code,
		     const int type_size,
		     const double * vdw_data,
		     double * ene_vdw,
		     double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double shX,shY,shZ;
  double jx,jy,jz;
  double d2;
  double d2inv;
  double d6inv;
  double d12inv;
  double shiftf;
  double switchf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double c2inv = 1.0/c2;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*d2inv;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;	  
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}


void 
pbcenergy_sh_rsswn_f(const int * first,
		     const int * second,
		     const int * boundarry,
		     const int natoms,
		     const float * coor,
		     const float * sh,
		     const float  c2,
		     const float  c_on2,
		     const float * charge,
		     const int * type_code,
		     const int type_size,
		     const float * vdw_data,
		     float * ene_vdw,
		     float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float shX,shY,shZ;
  float jx,jy,jz;
  float d2;
  float d2inv;
  float d6inv;
  float d12inv;
  float shiftf;
  float switchf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float c2inv = 1.0/c2;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	shiftf = one - d2*c2inv;
	shiftf = shiftf*shiftf;
	ecul = ecul + q1*q2*shiftf*d2inv;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  evdw  = evdw +
	    (nb12 - nb6)*switchf;	  
	}
	else {
	  evdw  = evdw + (nb12 - nb6);
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

//rswsn
void 
pbcenergy_rswsn_d(const int * first,
		  const int * second,
		  const int * boundarry,
		  const int natoms,
		  const double * coor,
		  const double  c2,
		  const double  c_on2,
		  const double * charge,
		  const int * type_code,
		  const int type_size,
		  const double * vdw_data,
		  double * ene_vdw,
		  double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double jx,jy,jz;
  double d2;
  double d6;
  double d2inv;
  double d6inv;
  double d12inv;
  double shiftf;
  double switchf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double c2inv = 1.0/c2;
  const double c6inv = c2inv*c2inv*c2inv;
  const double c12inv = c6inv*c6inv;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  ecul = ecul + q1*q2*switchf*d2inv;
	}
	else {
	  ecul = ecul + q1*q2*d2inv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

void 
pbcenergy_rswsn_f(const int * first,
		  const int * second,
		  const int * boundarry,
		  const int natoms,
		  const float * coor,
		  const float  c2,
		  const float  c_on2,
		  const float * charge,
		  const int * type_code,
		  const int type_size,
		  const float * vdw_data,
		  float * ene_vdw,
		  float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float jx,jy,jz;
  float d2;
  float d6;
  float d2inv;
  float d6inv;
  float d12inv;
  float shiftf;
  float switchf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float c2inv = 1.0/c2;
  const float c6inv = c2inv*c2inv*c2inv;
  const float c12inv = c6inv*c6inv;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  ecul = ecul + q1*q2*switchf*d2inv;
	}
	else {
	  ecul = ecul + q1*q2*d2inv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

void 
pbcenergy_sh_rswsn_d(const int * first,
		     const int * second,
		     const int * boundarry,
		     const int natoms,
		     const double * coor,
		     const double * sh,
		     const double  c2,
		     const double  c_on2,
		     const double * charge,
		     const int * type_code,
		     const int type_size,
		     const double * vdw_data,
		     double * ene_vdw,
		     double * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  double dx,dy,dz;
  double ix,iy,iz;
  double shX,shY,shZ;
  double jx,jy,jz;
  double d2;
  double d6;
  double d2inv;
  double d6inv;
  double d12inv;
  double shiftf;
  double switchf;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0, ecul=0;

  const double c2inv = 1.0/c2;
  const double c6inv = c2inv*c2inv*c2inv;
  const double c12inv = c6inv*c6inv;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  ecul = ecul + q1*q2*switchf*d2inv;
	}
	else {
	  ecul = ecul + q1*q2*d2inv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}

void 
pbcenergy_sh_rswsn_f(const int * first,
		     const int * second,
		     const int * boundarry,
		     const int natoms,
		     const float * coor,
		     const float * sh,
		     const float  c2,
		     const float  c_on2,
		     const float * charge,
		     const int * type_code,
		     const int type_size,
		     const float * vdw_data,
		     float * ene_vdw,
		     float * ene_cul){
  int f;
  int fpos;
  int s;
  int spos;
  float dx,dy,dz;
  float ix,iy,iz;
  float shX,shY,shZ;
  float jx,jy,jz;
  float d2;
  float d6;
  float d2inv;
  float d6inv;
  float d12inv;
  float shiftf;
  float switchf;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0, ecul=0;

  const float c2inv = 1.0/c2;
  const float c6inv = c2inv*c2inv*c2inv;
  const float c12inv = c6inv*c6inv;
  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;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    fpos = COOR_DIM*f;
    int b1 = boundarry[2*i];
    int b2 = boundarry[2*i+1];
    shX = sh[COOR_DIM*i];
    shY = sh[COOR_DIM*i+1];
    shZ = sh[COOR_DIM*i+2];
    ix = coor[fpos]+shX;
    iy = coor[fpos+1]+shY;
    iz = coor[fpos+2]+shZ;
    q1 = charge[f];
    tc1 = type_code[f];
    for(j= b1;j<b2;j++){
      s = second[j];
      spos = COOR_DIM*s;
      jx = coor[spos];
      jy = coor[spos+1];
      jz = coor[spos+2];
      dx = ix-jx;
      dy = iy-jy;
      dz = iz-jz;
      d2 = dx*dx+dy*dy+dz*dz;
      if(d2>c2) continue;
      {
	//compute energy
	q2 = charge[s];
	tc2 = type_code[s];
	data_pos = 4*(type_size*tc1+tc2);
	A = vdw_data[data_pos];
	B = vdw_data[data_pos+1];
	//now energy
	d2inv = one/d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb12 = A*(d12inv);
	nb6 = B*(d6inv);
	d6 = d2*d2*d2;
	shiftf = one - d6*c6inv;
	nb12 = A*(d12inv - c12inv*(one + shiftf+shiftf));
	nb6 = B*(d6inv -c6inv*(one + shiftf));
	evdw  = evdw +
	  nb12 - nb6;
	if(d2>c_on2){
	  switchf = c2diffinv*(c2-d2)*(c2-d2)*(c2+two*d2-three*c_on2);
	  ecul = ecul + q1*q2*switchf*d2inv;
	}
	else {
	  ecul = ecul + q1*q2*d2inv;
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
  *ene_cul = *ene_cul+ ecul;
}