/***************************************************************************
                          enefuncsoft.c  -  description
                             -------------------
    begin                : Fri Mar 18 12:25:31 2005
    copyright            : (C) 2002 by Cavalli Andrea
    email                : cavalli@bioc.unizh.ch
**************************************************************************/
/***************************************************************************
 *                                                                         *
 *   This program is free software; you can redistribute it and/or modify  *
 *   it under the terms of the GNU General Public License as published by  *
 *   the Free Software Foundation; either version 2 of the License, or     *
 *   (at your option) any later version.                                   *
 *                                                                         *
 ***************************************************************************/
#define COOR_DIM 4




//softswn
void 
energy_softswn_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  alpha,
		 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;
  double B;
  int tc1,tc2;
  int data_pos;
  double nb6;
  double nb12;
  double evdw =0;

  const double c2diffinv = 1.0/((c2-c_on2)*(c2-c_on2)*(c2-c_on2));
  const double one = 1.0;
  const double two = 2.0;
  const double three = 3.0;
  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];
    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
	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+alpha);
	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;	  
	}
	else {
	  evdw  = evdw + (nb12-nb6);
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
}



//softswn_f
void 
energy_softswn_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  alpha,
		 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;
  float B;
  int tc1,tc2;
  int data_pos;
  float nb6;
  float nb12;
  float evdw =0;

  const float c2diffinv = 1.0/((c2-c_on2)*(c2-c_on2)*(c2-c_on2));
  const float one = 1.0;
  const float two = 2.0;
  const float three = 3.0;
  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];
    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
	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+alpha);
	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;	  
	}
	else {
	  evdw  = evdw + (nb12-nb6);
	}
      }
    }
  }
      
  *ene_vdw = *ene_vdw + evdw;
}