// This file is part of the ESPResSo distribution (http://www.espresso.mpg.de).
// It is therefore subject to the ESPResSo license agreement which you accepted upon receiving the distribution
// and by which you are legally bound while utilizing this file in any form or way.
// There is NO WARRANTY, not even the implied warranty of MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
// You should have received a copy of that license along with this program;
// if not, refer to http://www.espresso.mpg.de/license.html where its current version can be found, or
// write to Max-Planck-Institute for Polymer Research, Theory Group, PO Box 3148, 55021 Mainz, Germany.
// Copyright (c) 2002-2006; all rights reserved unless otherwise stated.
#ifndef SUBT_GB_H
#define SUBT_GB_H

/** \file subt_gb.h
 *  Routines to calculate and subtract the Gay-Berne energy, torque and force 
 *  for a pair of particles.
 *  \ref forces.c
 *
 *  <b>Responsible:</b>
 *  <a href="mailto:address@hidden">Dmytro</a>
*/

#ifdef ROTATION

/// set the parameters for the subtract GB potential
MDINLINE int subt_gb_set_params(int bond_type, double k, double r)
{
  if(bond_type < 0)
    return TCL_ERROR;

  make_bond_type_exist(bond_type);

  bonded_ia_params[bond_type].p.subt_gb.k = k;
  bonded_ia_params[bond_type].p.subt_gb.r = r;
  bonded_ia_params[bond_type].type = BONDED_IA_SUBT_GB;  
  bonded_ia_params[bond_type].p.subt_gb.r2 = SQR(bonded_ia_params[bond_type].p.subt_gb.r);
  bonded_ia_params[bond_type].num = 1;

  mpi_bcast_ia_params(bond_type, -1); 

  return TCL_OK;
}

/// parse parameters for the subt_lj potential
MDINLINE int inter_parse_subt_gb(Tcl_Interp *interp, int bond_type, int argc, char **argv)
{
  double k, r;
  if (argc != 3) {
    Tcl_AppendResult(interp, "subt_gb needs 2 dummy parameters: "
		     "<k_subt_lj> <r_subt_lj>", (char *) NULL);
    return TCL_ERROR;
  }

  if ((! ARG_IS_D(1, k)) || (! ARG_IS_D(2, r))) {
    Tcl_AppendResult(interp, "subt_gb needs 2 dummy DOUBLE parameters: "
		     "<k_subt_lj> <r_subt_lj>", (char *) NULL);
    return TCL_ERROR;
  }

  CHECK_VALUE(subt_gb_set_params(bond_type, k, r), "bond type must be nonnegative");
}

MDINLINE int calc_subt_gb_pair_force(Particle *p1, Particle *p2, Bonded_ia_parameters *iaparams,
				double d[3], double force[3], double torque1[3], double torque2[3])

{  
   double Brack,BrackCut, Bra6,Bra6Cut, Bra12,Bra12Cut, Bra7,Bra7Cut, Bra13,Bra13Cut;
   double u1x, u1y, u1z,u2x, u2y, u2z, a, b, c;
   double chi2,chip2,X, Xcut;	  
   double E1,E2,E1E2, Sigma,soft;
   double apb,amb,apb2,amb2; 
   double dU_dr, dU_da, dU_db, dU_dc;  /*  all derivatives        */
   double denp,denm,denpp,denmp,den2,rep,repp,rem,remp,tmp,tmpp,den;
   double FikX,FikY,FikZ,dU_dc1,dU_dc2,dU_dc3;	        
/*  help for forces        */
//   double Gx,Gy,Gz;			
/*  help for torques       */

   IA_parameters *ia_params;
  double dist2 = sqrlen(d);
  double dist = sqrt(dist2);
  ia_params = get_ia_param(p1->p.type,p2->p.type);
 
  if(dist >= iaparams->p.subt_gb.r) return 1;


 if (dist < ia_params->GB_cut && ia_params->GB_eps != 0) {  

   u1x=p1->r.quatu[0];u1y=p1->r.quatu[1];u1z=p1->r.quatu[2];
   u2x=p2->r.quatu[0];u2y=p2->r.quatu[1];u2z=p2->r.quatu[2]; 
    
   a = (d[0]*u1x+ d[1]*u1y+d[2]*u1z)/dist;
   b = (d[0]*u2x+ d[1]*u2y+d[2]*u2z)/dist;
   c = u1x*u2x + u1y*u2y + u1z*u2z;	
   chi2   = ia_params->GB_chi1*ia_params->GB_chi1;
   chip2  = ia_params->GB_chi2*ia_params->GB_chi2;
   tmp = ia_params->GB_chi1*c; denp  = 1.0 + tmp ; denm  = 1.0 - tmp ;
   tmp = ia_params->GB_chi2*c; denpp = 1.0 + tmp ; denmp = 1.0 - tmp ;
   den2  = 1.0 - chi2*c*c;
   apb = a+b;
   amb = a-b;
   apb2=apb*apb;
   amb2=amb*amb;
  /******************************************************/
  /* Gay-Berne epsilon function */

   E1 = 1. / sqrt( den2 ) ;

  /* Gay-Berne epsilon prime function */

   E2 = 1.0 - 0.5 * ia_params->GB_chi2 * (apb2/denpp + amb2/denmp ) ;

  /* Gay-Berne well */

   if (ia_params->GB_mu == 2 && ia_params->GB_nu == 1) {
   E1E2 = ia_params->GB_eps*E1*E2*E2; 
   } else { 
   E1E2 = ia_params->GB_eps*pow(E1,ia_params->GB_nu)*pow(E2,ia_params->GB_mu);
   }
//  well = POWNU(eps) * POWMU(epsp) ;

  /* Gay-Berne sigma function */

   den = sqrt(1.0- 0.5* ia_params->GB_chi1 *(apb2/denp+amb2/denm )) ;
   Sigma = ia_params->GB_sig/den;
    soft=ia_params->GB_sig*1.00;
   X    = soft/(dist              - Sigma + soft);
   Xcut = soft/(ia_params->GB_cut - Sigma + soft);
     Bra6 = X*X*X;
     Bra6 = Bra6*Bra6;
     Bra7 = Bra6*X;
     Bra12 = Bra6*Bra6;
     Bra13 = Bra12*X;
     Brack = Bra6*(Bra6-1);
     if (ia_params->GB_cut > 0.0) 
     {Bra6Cut = Xcut*Xcut*Xcut;
      Bra6Cut = Bra6Cut*Bra6Cut;
      Bra7Cut = Bra6Cut*Xcut;
      Bra12Cut = Bra6Cut*Bra6Cut;
      Bra13Cut = Bra12Cut*Xcut;
      BrackCut = Bra6Cut*(Bra6Cut-1);
      Bra7  = Bra7-Bra7Cut;
      Bra13 = Bra13-Bra13Cut;
      Brack = Brack-BrackCut;}
     
   if (X < 1.5) { /* 1.25 corresponds to the interparticle penetration of 0.2 units of length.
		      If they are not that close, the GB forces and torques are calculated */
      /***************************************************/
      /* forces and torques between pair of molecules */

      rep  = (apb)/denp ;
      rem  = (amb)/denm ;
      repp = (apb)/denpp ;
      remp = (amb)/denmp ;

      /* derivative of potential energy with respect to molecular distance r(12) */
      
      tmp   = 3.0 * E1E2 * ( 2.0 * Bra13 - Bra7) ;
      dU_dr = - 2.0 * tmp /ia_params->GB_sig;
      tmp   *= ia_params->GB_chi1 / (den*den*den );

      /* derivative of potential energy 
	 with respect to cos(theta1) */

      tmpp = E1E2*ia_params->GB_mu*Brack*ia_params->GB_chi2/E2;
      dU_da = -tmpp*(repp+remp) + tmp*(rep+rem) ;

      /* derivative of potential energy with respect to cos(theta2) */

      dU_db = -tmpp*(repp-remp) + tmp*(rep-rem) ;

      /* derivative of potential energy with respect to cos(phi12) */

      dU_dc = E1E2*Brack*(ia_params->GB_nu*chi2*c/den2
        + 0.5*ia_params->GB_mu*chip2*(remp*remp-repp*repp)/E2)
        - 0.5*tmp*ia_params->GB_chi1*(rem*rem-rep*rep);

      /* rescale */
      dU_dr = 4*dU_dr/dist ;
      dU_da = 4*dU_da/dist ;
      dU_db = 4*dU_db/dist ;
      dU_dc = 4*dU_dc*1.;
      dU_dc1 = dU_dc*(u1y*u2z-u1z*u2y);
      dU_dc2 = dU_dc*(u1z*u2x-u1x*u2z);
      dU_dc3 = dU_dc*(u1x*u2y-u1y*u2x);
      a     /= dist ;
      b     /= dist ;

      /* force on molecule 1 due to molecule 2 */
      force[0] = (dU_dr*d[0]+dU_da*(u1x-a*d[0])+dU_db*(u2x-b*d[0]));
      force[1] = (dU_dr*d[1]+dU_da*(u1y-a*d[1])+dU_db*(u2y-b*d[1]));
      force[2] = (dU_dr*d[2]+dU_da*(u1z-a*d[2])+dU_db*(u2z-b*d[2]));
      
      /* torque on molecule 1 due to molecule 2 */
      torque1[0] = (dU_da*(u1y*d[2]-u1z*d[1])-dU_dc1);
      torque1[1] = (dU_da*(u1z*d[0]-u1x*d[2])-dU_dc2);
      torque1[2] = (dU_da*(u1x*d[1]-u1y*d[0])-dU_dc3);

      /* torque on molecule 2 due to molecule 1 */      
      torque2[0] = (dU_db*(u2y*d[2]-u2z*d[1])+dU_dc1);
      torque2[1] = (dU_db*(u2z*d[0]-u2x*d[2])+dU_dc2);
      torque2[2] = (dU_db*(u2x*d[1]-u2y*d[0])+dU_dc3);
            
   }
   else { /* the particles are too close to each other */
      /* derivative of potential energy with respect to molecular distance r(12) */      
      tmp   = 3.0 * E1E2 * X * ( 2.0 * Bra12 - Bra6 ) ;
      dU_dr = - 2.0 * tmp / ia_params->GB_sig;
     if (dU_dr > 0.) {dU_dr = 10/dist;} else {dU_dr = -10/dist;}
     FikX = dU_dr * d[0];
     FikY = dU_dr * d[1];
     FikZ = dU_dr * d[2];
     force[0] = FikX;
     force[1] = FikY;
     force[2] = FikZ;
 
 
  }
 }
// if (p1->p.identity==0 && p2->p.identity == 1) {printf("sub %f %f %f \n",torque1[0],torque1[1],torque1[2]);}
// if (p1->p.identity==1 && p2->p.identity == 0) {printf("sub %f %f %f \n",torque1[0],torque1[1],torque1[2]);}
// if (p1->p.identity==1 && p2->p.identity == 2) {printf("sub12 %f %f %f \n",FikX,FikY,FikZ);}
// if (p1->p.identity==2 && p2->p.identity == 1) {printf("sub21 %f %f %f \n",FikX,FikY,FikZ);}
   return 0;
 
}

MDINLINE double subt_gb_pair_energy(Particle *p1, Particle *p2, Bonded_ia_parameters *iaparams,
			       double d[3], double *_energy)
{double a,b,c, X, Xcut,
   Brack,BrackCut,
   u1x, u1y, u1z, u2x, u2y, u2z,	  
   E,E1,E2,soft, 
   Sigma,tmp,chi2,denp,denpp,denm,denmp,apb2,amb2,apb,amb;

   IA_parameters *ia_params;
   double dist2 = sqrlen(d);
   double dist = sqrt(dist2);
	
   if(dist >= iaparams->p.subt_gb.r) return 1;
   ia_params = get_ia_param(p1->p.type,p2->p.type);

   u1x=p1->r.quatu[0];u1y=p1->r.quatu[1];u1z=p1->r.quatu[2];
   u2x=p2->r.quatu[0];u2y=p2->r.quatu[1];u2z=p2->r.quatu[2]; 
    
   a = (d[0]*u1x+ d[1]*u1y+d[2]*u1z)/dist;
   b = (d[0]*u2x+ d[1]*u2y+d[2]*u2z)/dist;
   c = u1x*u2x + u1y*u2y + u1z*u2z;	
   chi2   = ia_params->GB_chi1*ia_params->GB_chi1;
   tmp = ia_params->GB_chi1*c; denp  = 1.0 + tmp ; denm  = 1.0 - tmp ;
   tmp = ia_params->GB_chi2*c; denpp = 1.0 + tmp ; denmp = 1.0 - tmp ;
   apb = a+b;    amb = a-b;
   apb2=apb*apb; amb2=amb*amb;
   E1 = 1. / sqrt(1.0 - chi2*c*c) ;
   E2 = 1.0 - 0.5 * ia_params->GB_chi2 * (apb2/denpp + amb2/denmp ) ;
   E = 4*ia_params->GB_eps*pow(E1,ia_params->GB_nu)*pow(E2,ia_params->GB_mu); 
   Sigma = ia_params->GB_sig/sqrt(1.0-0.5*ia_params->GB_chi1*(apb2/denp+amb2/denm));
   soft = ia_params->GB_sig*1.00;
   X    = soft/(dist              - Sigma + soft);
   Xcut = soft/(ia_params->GB_cut - Sigma + soft);

   Brack = X*X*X;
   Brack *= Brack;
   Brack = Brack*Brack-Brack;
     if (ia_params->GB_cut > 0.0) 
     {
      BrackCut = Xcut*Xcut*Xcut;
      BrackCut *= BrackCut;
      BrackCut = BrackCut*BrackCut-BrackCut;
      Brack = Brack-BrackCut;
     }

   *_energy = - E*Brack;
//   if (p1->p.identity == 0 && p2->p.identity== 1){printf("subtracting: %f \n",(-E*(Brack-BrackCut)) );}
//   if (p1->p.identity == 1 && p2->p.identity== 0){printf("subtracting: %f \n",(-E*(Brack-BrackCut)) );}
    return 0;
}

#endif
#endif