/***************************************************************************
                          sasaforcefunc.cpp  -  description
                             -------------------
    begin                : Fri Oct 15 14:02:17 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>
# define M_PI           3.14159265358979323846

void force1_rsssasa_d(const int * first,
		      const int * second,
		      const int * boundarry,
		      const int natoms,
		      const double * coor,
		      double * force,
		      const double  c2,
		      const double * charge,
		      const int * type_code,
		      const int type_size,
		      const double * vdw_data,
		      const int * sasa_ignore,
		      const double * sasa_vdw,
		      const double * sasa_surf,
		      const double * sasa_solv,
		      double * surf,
		      double  r_water){

  int f;
  int fpos;
  int s;
  int spos;
  int b1,b2;
  double fracX;
  double fracY;
  double fracZ;
  double dx,dy,dz;
  double ix,iy,iz;
  double jx,jy,jz;
  double tx,ty,tz;
  double d;
  double d2;
  double d6;
  double dinv;
  double d2inv;
  double d6inv;
  double d12inv;
  double sh;
  double frac;
  double nb6;
  double nb12;
  double nbelec;
  double A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  double si;
  double sj;
  double pi;
  double pj;
  double pij;
  double ri;
  double rj;
  double bi;
  double bj;
  double tmp;
  const double c2inv = 1.0/c2;
  const double c6inv = c2inv*c2inv*c2inv;
  const double c12inv = c6inv*c6inv;
  const double c18inv = c12inv*c6inv;
  const double one = 1.0;
  const double two = 2.0;
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    b1 = boundarry[2*i];
    b2 = boundarry[2*i+1];
    fracX = 0;
    fracY = 0;
    fracZ = 0;
    fpos = COOR_DIM*f;
    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){
	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];
	d2inv = one/d2;
	d6 = d2*d2*d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb6 = B*(d6inv-d6*c12inv); //mult by B
	nb12= A*(d12inv-d6*c18inv);
	frac = (nb12-nb6);
	sh = 1-d2*c2inv;
	nbelec = two*q1*q2*sh*(sh*d2inv+two*c2inv);
	frac = (frac + nbelec)*d2inv;
	tx = frac*dx;
	ty = frac*dy;
	tz = frac*dz;
	fracX = fracX - tx;
	fracY = fracY - ty;
	fracZ = fracZ - tz;
	force[spos] = force[spos] + tx;
	force[spos+1] = force[spos+1] + ty;
	force[spos+2] = force[spos+2] + tz;
	//begin sasa
	if(sasa_ignore[f]||
	   sasa_ignore[s]) continue;
	ri = sasa_vdw[f];
	rj = sasa_vdw[s];
	tmp = ri+rj+two*r_water;
	if((tmp*tmp)>= d2){
	  dinv = invsqrt(d2);
	  d = d2*dinv;
	  si = sasa_surf[f];
	  sj = sasa_surf[s];
	  pi = sasa_solv[f];
	  pj = sasa_solv[s];
	  pij = 0.3516;
	  bi = M_PI*(ri+r_water);
	  bj = M_PI*(rj+r_water);
	  tmp = tmp -d;
	  bi = bi*tmp;
	  bj = bj*tmp;
	  tmp = (rj-ri)*dinv;
	  bi = bi*(one+tmp);
	  bj = bj*(one-tmp);
	  surf[f] = surf[f]*(one-pi*pij*bi/si);
	  surf[s] = surf[s]*(one-pj*pij*bj/sj);
	}	
      }
    }
    force[fpos] = force[fpos] + fracX;
    force[fpos+1] = force[fpos+1] + fracY;
    force[fpos+2] = force[fpos+2] + fracZ;
  }
}


void force2_rsssasa_d(const int * first,
		      const int * second,
		      const int * boundarry,
		      const int natoms,
		      const double * coor,
		      double * force,
		      const int * sasa_ignore,
		      const double * sasa_vdw,
		      const double * sasa_surf,
		      const double * sasa_solv,
		      const double * sasa_sigma,
		      double * surf,
		      double  r_water){

  int f;
  int fpos;
  int s;
  int spos;
  int b1,b2;
  double fracX;
  double fracY;
  double fracZ;
  double frac;
  double dx,dy,dz;
  double ix,iy,iz;
  double jx,jy,jz;
  double tx,ty,tz;
  double d;
  double d2;
  double dinv;
  double d2inv;
  double si;
  double sj;
  double pi;
  double pj;
  double ri;
  double rj;
  double bi;
  double bj;
  double tmp;
  double tmp2;
  double bhat_i;
  double bhat_j;
  const double one = 1.0;
  const double two = 2.0;
  const double pij = 0.3516; 
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    if(sasa_ignore[f]) continue;
    b1 = boundarry[2*i];
    b2 = boundarry[2*i+1];
    fracX = 0;
    fracY = 0;
    fracZ = 0;
    fpos = COOR_DIM*f;
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    for(j= b1;j<b2;j++){
      s = second[j];
      if(sasa_ignore[s]) continue;
      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;
      ri = sasa_vdw[f];
      rj = sasa_vdw[s];
      tmp = ri+rj+two*r_water;
      if((tmp*tmp)>= d2){
	dinv = invsqrt(d2);
	d2inv = dinv*dinv;
	d = d2*dinv;
	si = sasa_surf[f];
	sj = sasa_surf[s];
	pi = sasa_solv[f];
	pj = sasa_solv[s];
	bi  =  M_PI*(ri+r_water); 
	tmp2 = tmp*(rj-ri)*d2inv;
	bhat_i = -bi; 
	bhat_i *= one+tmp2; 

	bj  =  M_PI*(rj+r_water); 
	bhat_j = -bj; 
	bhat_j *= one-tmp2;  

	tmp = tmp-d;
	tmp2 = (rj-ri)*dinv;
	bi *=  tmp; 
	bi *=  one+tmp2; 

	bj *=  tmp; 
	bj *=  one-tmp2; 	    
	bi = one/(si/(pi*pij)-bi)*bhat_i*dinv; 
	bj = one/(sj/(pj*pij)-bj)*bhat_j*dinv; 
	frac = sasa_sigma[f]*surf[f]*bi+sasa_sigma[s]*surf[s]*bj;
	tx = frac*dx;
	ty = frac*dy;
	tz = frac*dz;
	fracX = fracX - tx;
	fracY = fracY - ty;
	fracZ = fracZ - tz;
	force[spos] = force[spos] + tx;
	force[spos+1] = force[spos+1] + ty;
	force[spos+2] = force[spos+2] + tz;
      }
    }
    force[fpos] = force[fpos] + fracX;
    force[fpos+1] = force[fpos+1] + fracY;
    force[fpos+2] = force[fpos+2] + fracZ;
  }
}


void force1_rsssasa_f(const int * first,
		      const int * second,
		      const int * boundarry,
		      const int natoms,
		      const float * coor,
		      float * force,
		      const float  c2,
		      const float * charge,
		      const int * type_code,
		      const int type_size,
		      const float * vdw_data,
		      const int * sasa_ignore,
		      const float * sasa_vdw,
		      const float * sasa_surf,
		      const float * sasa_solv,
		      float * surf,
		      float  r_water){

  int f;
  int fpos;
  int s;
  int spos;
  int b1,b2;
  float fracX;
  float fracY;
  float fracZ;
  float dx,dy,dz;
  float ix,iy,iz;
  float jx,jy,jz;
  float tx,ty,tz;
  float d;
  float d2;
  float d6;
  float dinv;
  float d2inv;
  float d6inv;
  float d12inv;
  float sh;
  float frac;
  float nb6;
  float nb12;
  float nbelec;
  float A,B,q1,q2;
  int tc1,tc2;
  int data_pos;
  float si;
  float sj;
  float pi;
  float pj;
  float pij;
  float ri;
  float rj;
  float bi;
  float bj;
  float tmp;
  const float c2inv = 1.0/c2;
  const float c6inv = c2inv*c2inv*c2inv;
  const float c12inv = c6inv*c6inv;
  const float c18inv = c12inv*c6inv;
  const float one = 1.0;
  const float two = 2.0;
  int i,j; 
  for(i = 0;i<natoms;i++){
    f = first[i];
    b1 = boundarry[2*i];
    b2 = boundarry[2*i+1];
    fracX = 0;
    fracY = 0;
    fracZ = 0;
    fpos = COOR_DIM*f;
    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){
	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];
	d2inv = one/d2;
	d6 = d2*d2*d2;
	d6inv = d2inv*d2inv*d2inv;
	d12inv = d6inv*d6inv;
	nb6 = B*(d6inv-d6*c12inv); //mult by B
	nb12= A*(d12inv-d6*c18inv);
	frac = (nb12-nb6);
	sh = 1-d2*c2inv;
	nbelec = two*q1*q2*sh*(sh*d2inv+two*c2inv);
	frac = (frac + nbelec)*d2inv;
	tx = frac*dx;
	ty = frac*dy;
	tz = frac*dz;
	fracX = fracX - tx;
	fracY = fracY - ty;
	fracZ = fracZ - tz;
	force[spos] = force[spos] + tx;
	force[spos+1] = force[spos+1] + ty;
	force[spos+2] = force[spos+2] + tz;
	//begin sasa
	if(sasa_ignore[f]||
	   sasa_ignore[s]) continue;
	ri = sasa_vdw[f];
	rj = sasa_vdw[s];
	tmp = ri+rj+two*r_water;
	if((tmp*tmp)>= d2){
	  dinv = invsqrt(d2);
	  d = d2*dinv;
	  si = sasa_surf[f];
	  sj = sasa_surf[s];
	  pi = sasa_solv[f];
	  pj = sasa_solv[s];
	  pij = 0.3516;
	  bi = M_PI*(ri+r_water);
	  bj = M_PI*(rj+r_water);
	  tmp = tmp -d;
	  bi = bi*tmp;
	  bj = bj*tmp;
	  tmp = (rj-ri)*dinv;
	  bi = bi*(one+tmp);
	  bj = bj*(one-tmp);
	  surf[f] = surf[f]*(one-pi*pij*bi/si);
	  surf[s] = surf[s]*(one-pj*pij*bj/sj);
	}	
      }
    }
    force[fpos] = force[fpos] + fracX;
    force[fpos+1] = force[fpos+1] + fracY;
    force[fpos+2] = force[fpos+2] + fracZ;
  }
}


void force2_rsssasa_f(const int * first,
		      const int * second,
		      const int * boundarry,
		      const int natoms,
		      const float * coor,
		      float * force,
		      const int * sasa_ignore,
		      const float * sasa_vdw,
		      const float * sasa_surf,
		      const float * sasa_solv,
		      const float * sasa_sigma,
		      float * surf,
		      float  r_water){

  int f;
  int fpos;
  int s;
  int spos;
  int b1,b2;
  float fracX;
  float fracY;
  float fracZ;
  float frac;
  float dx,dy,dz;
  float ix,iy,iz;
  float jx,jy,jz;
  float tx,ty,tz;
  float d;
  float d2;
  float dinv;
  float d2inv;
  float si;
  float sj;
  float pi;
  float pj;
  float ri;
  float rj;
  float bi;
  float bj;
  float tmp;
  float tmp2;
  float bhat_i;
  float bhat_j;
  const float one = 1.0;
  const float two = 2.0;
  const float pij = 0.3516; 
  int i,j;
  for(i = 0;i<natoms;i++){
    f = first[i];
    if(sasa_ignore[f]) continue;
    b1 = boundarry[2*i];
    b2 = boundarry[2*i+1];
    fracX = 0;
    fracY = 0;
    fracZ = 0;
    fpos = COOR_DIM*f;
    ix = coor[fpos];
    iy = coor[fpos+1];
    iz = coor[fpos+2];
    for(j= b1;j<b2;j++){
      s = second[j];
      if(sasa_ignore[s]) continue;
      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;
      ri = sasa_vdw[f];
      rj = sasa_vdw[s];
      tmp = ri+rj+two*r_water;
      if((tmp*tmp)>= d2){
	dinv = invsqrt(d2);
	d2inv = dinv*dinv;
	d = d2*dinv;
	si = sasa_surf[f];
	sj = sasa_surf[s];
	pi = sasa_solv[f];
	pj = sasa_solv[s];
	bi  =  M_PI*(ri+r_water); 
	tmp2 = tmp*(rj-ri)*d2inv;
	bhat_i = -bi; 
	bhat_i *= one+tmp2; 

	bj  =  M_PI*(rj+r_water); 
	bhat_j = -bj; 
	bhat_j *= one-tmp2;  

	tmp = tmp-d;
	tmp2 = (rj-ri)*dinv;
	bi *=  tmp; 
	bi *=  one+tmp2; 

	bj *=  tmp; 
	bj *=  one-tmp2; 	    
	bi = one/(si/(pi*pij)-bi)*bhat_i*dinv; 
	bj = one/(sj/(pj*pij)-bj)*bhat_j*dinv; 
	frac = sasa_sigma[f]*surf[f]*bi+sasa_sigma[s]*surf[s]*bj;
	tx = frac*dx;
	ty = frac*dy;
	tz = frac*dz;
	fracX = fracX - tx;
	fracY = fracY - ty;
	fracZ = fracZ - tz;
	force[spos] = force[spos] + tx;
	force[spos+1] = force[spos+1] + ty;
	force[spos+2] = force[spos+2] + tz;
      }
    }
    force[fpos] = force[fpos] + fracX;
    force[fpos+1] = force[fpos+1] + fracY;
    force[fpos+2] = force[fpos+2] + fracZ;
  }
}