__global__ void solvation (int atoms_r, int atoms_l, int nlig, float *d_rec_x, float *d_rec_y, float *d_rec_z, float *d_lig_x, float *d_lig_y, float *d_lig_z, int* d_rectype, int* d_ligtype, float *d_ql ,float *d_qr, float *d_energy, struct autodock_param_t *d_a_params, int nconformations) {

  float dist, temp_desolv = 0, miatomo[3], e_desolv;
  int ind1, ind2;
  int total;

  float difx, dify, difz, solv_asp_1, solv_asp_2, solv_vol_1, solv_vol_2, solv_qasp_1, solv_qasp_2;
  float mod2x, mod2y, mod2z;

  unsigned int i = blockIdx.x + blockDim.x + threadIdx.x;
  int conformacion = i / nlig;
  int atomo = i % nlig;
  int tid = threadIdx.x;

  e_desolv = 0;
  ind1 = d_ligtype[atomo];
  miatomo[0] = *(d_lig_x + conformacion*nlig + atomo);
  miatomo[1] = *(d_lig_y + conformacion*nlig + atomo);
  miatomo[2] = *(d_lig_z + conformacion*nlig + atomo);
  solv_asp_1 = d_a_params[ind1].asp;
  solv_vol_1 = d_a_params[ind1].vol;

  for(int j = 0; j < atoms_l; j++){

    ind2 = d_rectype[j];
    solv_asp_2 = d_a_params[ind2].asp;
    solv_vol_2 = d[atomo]ms[ind2].vol;
    difx = (d_rec_x[j]) - miatomo[0];
    dify = (d_rec_y[j]) - miatomo[1];
    difz = (d_rec_z[j]) - miatomo[2];
    mod2x = difx * difx;
    mod2y = dify * dify;
    mod2z = difz * difz;
    difx = mod2x + mod2y + mod2z; {
    e_desolv = ((solv_asp_1 * solv_vol_2) + (QASP * fabs(d_ql[atomo]) * solv_vol_2) + (solv_asp_2 * solv_vol_1) + (QASP * fabs(d_qr[j]) * solv_vol_1) * exp(-difx/(2*G_D_2)));
    printf("%f\n",e_desolv);
   }

   d_energy[i] = e_desolv;
}


void forces_GPU_AU (int atoms_r, int atoms_l, int nlig, float *rec_x, float *rec_y, float *rec_z, float *lig_x, float *lig_y, float *lig_z, int* rectype, int* ligtype, float *ql ,float *qr, float *energy, struct autodock_param_t *a_params, int nconformations){

  float *d_rec_x, *d_rec_y, *d_rec_z, *d_lig_x, *d_lig_y, *d_lig_z, *d_ql, *d_qr, *d_energy, aux=0;
  int *d_rectype, *d_ligtype,i;
  struct autodock_param_t *d_a_params;

  cudaSetDevice(0);
  cudaMalloc((void**)&d_rec_x, sizeof(rec_x));
  cudaMalloc((void**)&d_rec_y, sizeof(rec_y));
  cudaMalloc((void**)&d_rec_z, sizeof(rec_z));
  cudaMalloc((void**)&d_lig_x, sizeof(lig_x));
  cudaMalloc((void**)&d_lig_y, sizeof(lig_y));
  cudaMalloc((void**)&d_lig_z, sizeof(lig_z));
  cudaMalloc((void**)&d_ql, sizeof(ql));
  cudaMalloc((void**)&d_qr, sizeof(qr));
  cudaMalloc((void**)&d_energy, sizeof(energy));
  cudaMalloc((void**)&d_rectype, sizeof(rectype));
  cudaMalloc((void**)&d_ligtype, sizeof(ligtype));
  cudaMalloc((void**)&d_a_params, sizeof(a_params));

  cudaMemcpy(d_rec_x,rec_x,sizeof(rec_x),cudaMemcpyHostToDevice);
  cudaMemcpy(d_rec_y,rec_y,sizeof(rec_y),cudaMemcpyHostToDevice);
  cudaMemcpy(d_rec_z,rec_z,sizeof(rec_z),cudaMemcpyHostToDevice);
  cudaMemcpy(d_lig_x,lig_x,sizeof(lig_x),cudaMemcpyHostToDevice);
  cudaMemcpy(d_lig_y,lig_y,sizeof(lig_y),cudaMemcpyHostToDevice);
  cudaMemcpy(d_lig_z,lig_z,sizeof(lig_z),cudaMemcpyHostToDevice);
  cudaMemcpy(d_ql,ql,sizeof(ql),cudaMemcpyHostToDevice);
  cudaMemcpy(d_qr,qr,sizeof(qr),cudaMemcpyHostToDevice);
  cudaMemcpy(d_energy,energy,sizeof(energy),cudaMemcpyHostToDevice);
  cudaMemcpy(d_rectype,rectype,sizeof(rectype),cudaMemcpyHostToDevice);
  cudaMemcpy(d_ligtype,ligtype,sizeof(ligtype),cudaMemcpyHostToDevice);
  cudaMemcpy(d_a_params,a_params,sizeof(a_params),cudaMemcpyHostToDevice);

  //Definir numero de hilos y bloques
  dim3 thread (tam_block);
  dim3 block ((nconformations*nlig/tam_block)+1);

  //Anadir parametros
  solvation <<< block,thread>>> (atoms_r,atoms_l, nlig, d_rec_x,d_rec_y,d_rec_z,d_lig_x,d_lig_y,d_lig_z,d_rectype,d_ligtype,d_ql,d_qr,d_energy,d_a_params, nconformations);

  cudaMemcpy(energy,d_energy,sizeof(energy),cudaMemcpyDeviceToHost);
  //for(i = 0; i < (nconformations*nlig) ; i++){
    // aux += energy[i];
  //}
  //printf("aux = %f",aux);
}