# include "wpmd_split.h"
//# include "erf.h"


void AWPMD_split::resize(int flag){
  for(int s=0;s<2;s++){
    wf_norm[s].resize(ne[s]);
    if(flag&(0x8|0x4)){ //electron forces needed
      wf_norm_der[s].resize(nvar[s]);
      ovl_der[s].resize(nvar[s]);
      E_der[s].resize(nvar[s]);

      if(flag&(0x8|0x4) || norm_needed){ //electron forces or norm matrix needed 
        if(approx==HARTREE){ // L and Norm are needed in block form
          Lh[s].resize(nvar[s]);
          if(norm_needed){
            Normh[s].resize(ne[s]);
            for(int i=0;i<ne[s];i++)
              Normh[s][i].init(10*nspl[s][i]);
          }
        }
        else if(norm_needed){
          Norm[s].init(nvar[s]);
        }
      }
    }
  }  
  AWPMD::resize(flag);
}


int AWPMD_split::add_split(Vector_3 &x, Vector_3 &v, double &w, double &pw, Vector_2 &c, double mass, double q, int tag){
  if(!spl_add){
    nspl[s_add].push_back(1);
    ne[s_add]++;
  }
  else{
    nspl[s_add][ne[s_add]-1]++; // incrementing the WP number for the last electron
  }
  spl_add++;
  nwp[s_add]++;
  nvar[s_add]+=10;
  wp[s_add].push_back(create_wp(x,v,w,pw,mass));
  split_c[s_add].push_back(c);
  qe[s_add].push_back(q);
  
  if(tag==0)
    tag=spl_add;
  partition1[s_add].push_back(tag);
  
  return spl_add-1;
}

int AWPMD_split::set_electrons(int s, int nel, Vector_3P x, Vector_3P v, double* w, double* pw, Vector_2 *c, int *splits, double mass, double *q)
{
  if(s < 0 || s > 1)
    return LOGERR(-1,fmt("AWPMD_split.set_electrons: invaid spin setting (%d)!",s),LINFO);

  // calculating the total n
  nvar[s]=0;
  int n=0;
  for(int i=0;i<nel;i++){
    n+=splits[i];
    nvar[s]+=10*splits[i]; // number of dynamic variables per wp: x[3],p[3],w,pw,c_re,c_im
  }
  nwp[s]=n;

  norm_matrix_state[s] = NORM_UNDEFINED;
  ne[s]=nel;
  wp[s].resize(n);
  
  split_c[s].resize(n);
  split_c[s].assign(c,c+n);
  
  
  nspl[s].resize(nel);
  nspl[s].assign(splits,splits+nel);

  partition1[s].clear();
  for(int i=0;i<n;i++){
    
    /*if(constraint==FIX){
      w[i]=w0;
      pw[i]=0.;
    }
    
    double rw;
    if(Lextra>0){ // width PBC, keeping the width are within [0,Lextra]
      w[i]=fmod(w[i],Lextra);
      if(w[i]<0) w[i]+=Lextra;
      rw=w[i]; // WP width for energy evaluation is within [0, L/2]
      if(rw > Lextra/2) rw = Lextra - rw;
    }
    else
      rw=w[i];

    wp[s][i].init(rw,x[i],v[i]*m_electron/h_plank,pw[i]/h_plank);*/
    wp[s][i]=create_wp(x[i],v[i],w[i],pw[i],mass);
    //printf("%15d %15g\n",i,rw);
    // assign default partition
    partition1[s].push_back(i+1);
  }

  // assign electronic charge
  if(q)
    qe[s].assign(q,q+nwp[s]);
  else
    qe[s].assign(nwp[s],-1);

  

  return 1;
}



void AWPMD_split::eterm_deriv(int ic1,int s1,int c1, int j1,int ic2,int s2, int c2, int k2,cdouble pref,
                              const OverlapDeriv &o,cdouble v,cdouble dv_aj_conj,
                              cdouble dv_ak,cVector_3 dv_bj_conj, cVector_3 dv_bk){
  cdouble cj(split_c[s1][ic1+j1][0],split_c[s1][ic1+j1][1]);
  cdouble ck(split_c[s2][ic2+k2][0],split_c[s2][ic2+k2][1]);
  int indw1=8*ic1, indw2=8*ic2;
  int indn1=(nvar[s1]/10)*8+2*ic1, indn2=(nvar[s2]/10)*8+2*ic2;
  cdouble part_jk=conj(cj)*ck; 

  int M= 1; //(j1==k2 ? 1 : 2);

  // over a_k_re
  E_der[s2][indw2+8*k2]+=   M*real(pref*part_jk*(o.da2_re()*v+o.I0*dv_ak));           
  // over a_k_im
  E_der[s2][indw2+8*k2+1]+=   M*real(pref*part_jk*(o.da2_im()*v+i_unit*o.I0*dv_ak));    
  // over a_j_re
  E_der[s1][indw1+8*j1]+=   M*real(pref*part_jk*(o.da1_re()*v+o.I0*dv_aj_conj));           
  // over a_j_im
  E_der[s1][indw1+8*j1+1]+= M*real(pref*part_jk*(o.da1_im()*v-i_unit*o.I0*dv_aj_conj));           

  for(int i=0;i<3;i++){
    // over b_k_re
    E_der[s2][indw2+8*k2+2+2*i]+= M*real(pref*part_jk*(o.db2_re(i)*v+o.I0*dv_bk[i]));           
    // over b_k_im
    E_der[s2][indw2+8*k2+2+2*i+1]+= M*real(pref*part_jk*(o.db2_im(i)*v+i_unit*o.I0*dv_bk[i]));      
    // over b_j_re
    E_der[s1][indw1+8*j1+2+2*i]+= M*real(pref*part_jk*(o.db1_re(i)*v+o.I0*dv_bj_conj[i]));                      
    // over b_j_im
    E_der[s1][indw1+8*j1+2+2*i+1]+= M*real(pref*part_jk*(o.db1_im(i)*v-i_unit*o.I0*dv_bj_conj[i]));             
  }


  // over ck_re
  E_der[s2][indn2+2*k2]+=M*real(pref*conj(cj)*o.I0*v);  
  // over ck_im
  E_der[s2][indn2+2*k2+1]+=M*real(pref*i_unit*conj(cj)*o.I0*v);
  // over cj_re
  E_der[s1][indn1+2*j1]+=M*real(pref*ck*o.I0*v);  
  // over cj_im
  E_der[s1][indn1+2*j1+1]+=M*real(-pref*i_unit*ck*o.I0*v);  
  
  double t= -M*real(pref*part_jk*o.I0*v);
  // nonlocal terms: TODO: make a separate global loop for summation of nonlocal terms
  for(int j=0;j<nspl[s1][c1];j++){
    for(int i=0;i<8;i++)
      E_der[s1][indw1+8*j+i]+=t*wf_norm_der[s1][indw1+8*j+i];
    E_der[s1][indn1+2*j  ]+=t*wf_norm_der[s1][indn1+2*j  ];
    E_der[s1][indn1+2*j+1]+=t*wf_norm_der[s1][indn1+2*j+1];
  }
  for(int k=0;k<nspl[s2][c2];k++){
    for(int i=0;i<8;i++)
      E_der[s2][indw2+8*k+i]+=t*wf_norm_der[s2][indw2+8*k+i];
    E_der[s2][indn2+2*k]+=   t*wf_norm_der[s2][indn2+2*k  ];
    E_der[s2][indn2+2*k+1]+= t*wf_norm_der[s2][indn2+2*k+1];
  }

}


void AWPMD_split::calc_norms(int flag){
  if(flag&0x4){ // electron forces requested
    for(int s1=0;s1<2;s1++){ // clearing norm derivatives
      for(int i=0;i<nvar[s1];i++){
        wf_norm_der[s1][i]=0;
        E_der[s1][i]=0;
        ovl_der[s1][i]=0;
      }
    }
  }
  // calculating block norms and derivatives
  for(int s1=0;s1<2;s1++){
    int ic1=0; // starting index of the wp for current electron
    int indw1=0; // starting index of the electron wp coordinates 
    int indn1=(nvar[s1]/10)*8; // starting index of the electron norm coordinates 
    
    for(int c1=0;c1<ne[s1];c1++){
      
      // calculating the block norm
      wf_norm[s1][c1]=0.;
      for(int j1=0;j1<nspl[s1][c1];j1++){
        double cj_re=split_c[s1][ic1+j1][0];
        double cj_im=split_c[s1][ic1+j1][1];
        cdouble ccj=cdouble(cj_re,-cj_im);
        double part_jj=norm(ccj);
        wf_norm[s1][c1]+=part_jj;
        WavePacket& wj = wp[s1][ic1+j1];
        OverlapDeriv o;
        if(flag&(0x8|0x4)){ //electron forces needed
          wf_norm_der[s1][indn1+2*j1]+=2*cj_re;  // over cj_re
          wf_norm_der[s1][indn1+2*j1+1]+=2*cj_im; // over cj_im
          o.set1(wj);// conjugate: mu -> -mu, v -> -v !!!
        }
        
       
        for(int k1=j1+1;k1<nspl[s1][c1];k1++){
          double ck_re=split_c[s1][ic1+k1][0];
          double ck_im=split_c[s1][ic1+k1][1];
          cdouble ck=cdouble(ck_re,ck_im);

          WavePacket wk=wp[s1][ic1+k1];
          if(pbc) 
            move_to_image(wj,wk);
          WavePacket wjk=conj(wj)*wk;
          cdouble I0=wjk.integral();
          
          cdouble part_jk=ccj*ck;
          wf_norm[s1][c1]+=2*real(part_jk*I0);

          
          if(flag&(0x8|0x4)){ //electron forces needed
            o.set2(wk,&I0);
            
            
            wf_norm_der[s1][indw1+8*k1]+=   2*real(part_jk*o.da2_re());           // over a_k_re
            wf_norm_der[s1][indw1+8*k1+1]+=   2*real(part_jk*o.da2_im());    // over a_k_im

            wf_norm_der[s1][indw1+8*j1]+=   2*real(part_jk*o.da1_re());             // over a_j_re
            wf_norm_der[s1][indw1+8*j1+1]+= 2*real(part_jk*o.da1_im());    // over a_j_im

            for(int i=0;i<3;i++){
              wf_norm_der[s1][indw1+8*k1+2+2*i]+= 2*real(part_jk*o.db2_re(i));           // over b_k_re
              wf_norm_der[s1][indw1+8*k1+2+2*i+1]+= 2*real(part_jk*o.db2_im(i));  // over b_k_im

              wf_norm_der[s1][indw1+8*j1+2+2*i]+= 2*real(part_jk*o.db1_re(i));           // over b_j_re
              wf_norm_der[s1][indw1+8*j1+2+2*i+1]+= 2*real(part_jk*o.db1_im(i));  // over b_j_im
            }
     
            wf_norm_der[s1][indn1+2*j1]+=2*real(ck*I0);  // over cj_re
            wf_norm_der[s1][indn1+2*j1+1]+=2*real(-i_unit*ck*I0);  // over cj_im
            wf_norm_der[s1][indn1+2*k1]+=2*real(ccj*I0); // over ck_re
            wf_norm_der[s1][indn1+2*k1+1]+=2*real(i_unit*ccj*I0); // over ck_im

            
            // overlap derivatives (for norm matrix)
            ovl_der[s1][indw1+8*k1]+=   ccj*o.da2_re();           // over a_k_re
            ovl_der[s1][indw1+8*k1+1]+=   ccj*o.da2_im();    // over a_k_im

            ovl_der[s1][indw1+8*j1]+=   ck*o.da1_re();             // over a_j_re
            ovl_der[s1][indw1+8*j1+1]+= ck*o.da1_im();    // over a_j_im

            for(int i=0;i<3;i++){
              ovl_der[s1][indw1+8*k1+2+2*i]+= ccj*o.db2_re(i);           // over b_k_re
              ovl_der[s1][indw1+8*k1+2+2*i+1]+= ccj*o.db2_im(i);  // over b_k_im

              ovl_der[s1][indw1+8*j1+2+2*i]+= ck*o.db1_re(i);           // over b_j_re
              ovl_der[s1][indw1+8*j1+2+2*i+1]+= ck*o.db1_im(i);  // over b_j_im
            }
     
            ovl_der[s1][indn1+2*j1]+=ck*I0;  // over cj_re
            ovl_der[s1][indn1+2*j1+1]+=-i_unit*ck*I0;  // over cj_im
            ovl_der[s1][indn1+2*k1]+=ccj*I0; // over ck_re
            ovl_der[s1][indn1+2*k1+1]+=i_unit*ccj*I0; // over ck_im
          }
        } // k1
      }// j1
      // normalizing the norm derivative 
      for(int j1=0;j1<nspl[s1][c1];j1++){
        for(int i=0;i<8;i++) // wp parameters
          wf_norm_der[s1][indw1+8*j1+i]/=2*wf_norm[s1][c1];
        // c
        wf_norm_der[s1][indn1+2*j1]/=2*wf_norm[s1][c1];  
        wf_norm_der[s1][indn1+2*j1+1]/=2*wf_norm[s1][c1];
      }
      //wf_norm[s1][c1]=1; 
      //wf_norm[s1][c1]=sqrt(wf_norm[s1][c1]);
      ic1+=nspl[s1][c1];
      indw1+=8*nspl[s1][c1]; // 8 variables in each wavepacket
      indn1+=2*nspl[s1][c1]; // 2 variables in each wp norm
    }// c1   
  }// s1
}

void AWPMD_split::clear_forces(int flag,Vector_3P fi, Vector_3P fe_x, 
                               Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c){
  if(flag&0x1){
    for(int i=0;i<ni;i++)
      fi[i]=Vector_3(0.);
  }
  if(flag&0x4 && !(flag&0x10)){ // electron forces requested in physical representation
    int iv1=0;
    for(int s1=0;s1<2;s1++){ // clearing forces
      for(int c1=0;c1<ne[s1];c1++){
        for(int j1=0;j1<nspl[s1][c1];j1++){
          fe_x[iv1+j1]=Vector_3(0,0,0);
          fe_p[iv1+j1]=Vector_3(0,0,0);
          fe_w[iv1+j1]=0;
          fe_pw[iv1+j1]=0;
          fe_c[iv1+j1]=Vector_2(0,0);
        }
        iv1+=nspl[s1][c1];
      }
    }
  }
}

void AWPMD_split::get_el_forces(int flag, Vector_3P fe_x, 
                               Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c){
  if(flag&0x4) //need to replace the forces
    clear_forces(0x4,NULL,fe_x,fe_p,fe_w,fe_pw,fe_c);
  
  // recalculating derivatives
  if(flag&(0x8|0x4)){ //electron forces needed
    int iv1=0;
    for(int s1=0;s1<2;s1++){
      int ic1=0; // starting index of the wp for current electron
      for(int c1=0;c1<ne[s1];c1++){
        int indw1=8*ic1;
        int indn1=(nvar[s1]/10)*8+2*ic1;
        for(int k1=0;k1<nspl[s1][c1];k1++){
          WavePacket wk=wp[s1][ic1+k1];  
          /*double w=wk.get_width();
          Vector_3 r=wk.get_r();
          double t=3/(2*w*w*w);
          fe_w[ic1+k1]+= t*E_der[s1][indw1+8*k1]+imag(wk.a)*E_der[s1][indw1+8*k1+1]/w;
          fe_pw[ic1+k1]+=E_der[s1][indw1+8*k1+1]/(2*w*h_plank);
          for(int i=0;i<3;i++){
            fe_x[ic1+k1][i]+= -2*real(wk.a)*E_der[s1][indw1+8*k1+2+2*i]-2*imag(wk.a)*E_der[s1][indw1+8*k1+2+2*i+1];
            fe_p[ic1+k1][i]+= (-E_der[s1][indw1+8*k1+2+2*i+1])*(m_electron/h_plank); //*(h_plank/m_electron);
            fe_pw[ic1+k1]+=(r[i]*E_der[s1][indw1+8*k1+2+2*i+1]/w)/h_plank;  
            fe_w[ic1+k1]+=2*r[i]*(t*E_der[s1][indw1+8*k1+2+2*i]+imag(wk.a)*E_der[s1][indw1+8*k1+2+2*i+1]/w);
          }*/
          wk.int2phys_der< minus >(E_der[s1].begin()+indw1,(double *)&fe_x[iv1+k1],(double *)&fe_p[iv1+k1],&fe_w[iv1+k1],&fe_pw[iv1+k1], 1./one_h);
          fe_c[iv1+k1]+=-Vector_2(E_der[s1][indn1+2*k1],E_der[s1][indn1+2*k1+1]);
        }// k1
        ic1+=nspl[s1][c1]; // incrementing block1 wp address
        iv1+=nspl[s1][c1]; // incrementing global variable address
      }
    }  
  }
}

//e same as interaction, but using Hartee factorization (no antisymmetrization)
int  AWPMD_split::interaction_hartree(int flag, Vector_3P fi, Vector_3P fe_x, 
                                      Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c){
  
  // resize arrays if needed
  enum APPROX tmp=HARTREE;
  swap(tmp,approx); // do not neeed large matrices
  resize(flag);
  swap(tmp,approx);

  // clearing forces
  clear_forces(flag,fi,fe_x,fe_p,fe_w,fe_pw,fe_c);
  // calculate block norms and (optionally) derivatives
  calc_norms(flag);  

  Eee = Ew = 0.;
  for(int s1=0;s1<2;s1++){
    Ee[s1]=0.;
    Eei[s1]=0.;
    int ic1=0; // starting index of the wp for current electron
    
    for(int c1=0;c1<ne[s1];c1++){
      // calculating single-electron quantities within block
      double Ee1=0., Ew1=0., Eei1=0.;
      double pref=-h2_me/(2*wf_norm[s1][c1]); // ekin
      double pref_ei=coul_pref/wf_norm[s1][c1];
      
      for(int j1=0;j1<nspl[s1][c1];j1++){
        cdouble cj(split_c[s1][ic1+j1][0],split_c[s1][ic1+j1][1]);
        WavePacket wj=wp[s1][ic1+j1];  

        OverlapDeriv o;
        if(flag&(0x8|0x4)) //electron forces needed
          o.set1(wj);
          
        for(int k1=j1;k1<nspl[s1][c1];k1++){
          int M1a=(j1==k1 ? 1: 2);
          double M1e, M1f;
          _mytie(M1e,M1f)=check_part1(s1,ic1+j1,ic1+k1)*M1a;
          
          cdouble ck(split_c[s1][ic1+k1][0],split_c[s1][ic1+k1][1]);
          
          // electrons kinetic energy
          WavePacket wk=wp[s1][ic1+k1];
          if(pbc)
            move_to_image(wj,wk);
          
          WavePacket wjk=conj(wj)*wk;
          cdouble I0 = wjk.integral();
          cdouble part_jk=conj(cj)*ck;

          if(M1e){
            cVector_3 v1=conj(wj.b)*wk.a-wk.b*conj(wj.a);
            cdouble v=(v1*v1)/wjk.a;
            v-=6.*wk.a*conj(wj.a);
            v/=wjk.a;
            
            //v=1.; 
           
            // kinetic energy contribution
            Ee[s1] += M1e*real(part_jk*I0*v)*pref;  // Ejk+Ekj=Ejk+conj(Ejk), j!=k or Ejj
            


            if(flag&(0x8|0x4)){ //electron forces needed
              cVector_3 tv=wk.b*conj(wj.a)-conj(wj.b)*wk.a;
              cdouble ajk2=wjk.a*wjk.a;
              cdouble ajk3=ajk2*wjk.a;
              cdouble dv_aj_conj=     -2*wk.a*(3*wjk.a*wk.a-tv*wjk.b)/ajk3;
              cdouble dv_ak=-2*conj(wj.a)*(3*wjk.a*conj(wj.a)+tv*wjk.b)/ajk3;
              cVector_3 dv_bj_conj=-2*tv*wk.a/ajk2;
              cVector_3 dv_bk=2*tv*conj(wj.a)/ajk2;

              /*cdouble dv_aj_conj=0.;
              cdouble dv_ak=0.;
              cVector_3 dv_bj_conj;
              cVector_3 dv_bk;*/
            
                 
            
              o.set2(wk,&I0);
              // calculate full derivative of the term pref*conj(cj)*ck*Ijk*vjk/sqrt(nrm(s1)*nrm(s2))
              // denominator must be included in pref
              eterm_deriv(ic1,s1,c1,j1,ic1,s1,c1,k1,M1f*pref,o,v,dv_aj_conj,dv_ak,dv_bj_conj,dv_bk);
            }
          }// M1e

          cVector_3 djk=wjk.b/(2.*wjk.a);
          // e-i energy
          //cdouble sum(0.,0.);
          for(int i=0;i<ni;i++){  // ions loop 
            double M1ie, M1if;
            _mytie(M1ie,M1if)=check_part1ei(s1,ic1+j1,ic1+k1,i)*M1a;
            if(!M1ie)
              continue;

            cVector_3 gjki=djk - cVector_3(xi[i]);
            
            if(pbc) // correcting the real part (distance) according to PBC
              gjki=rcell1(gjki,cell,pbc);
              //-Igor- gkli=cVector_3(real(gkli).rcell1(cell,pbc),imag(gkli));
            
            cdouble ngjki=gjki.norm();
            cdouble sqajk=sqrt(wjk.a);
            //cdouble ttt = cerf_div(ngjki,c);
            cdouble v=cerf_div(ngjki,sqajk);
            double pref_eiq=pref_ei*qi[i]*(qe[s1][ic1+j1]+qe[s1][ic1+k1])*0.5;
            cdouble dE=pref_eiq*I0*v*part_jk;

            // ions energy contribution
            Eei[s1] += M1ie*real(dE);
            //sum+=dE;

            cdouble t1, t2;
            bool zero_force=(fabs(real(ngjki))+fabs(imag(ngjki))<1e-10);
            if(flag&(0x8|0x4|0x3) ){ //some derivatives needed
              cdouble arg=ngjki*sqajk;
              t1=two_over_sqr_pi*exp(-arg*arg);
              t2=t1/(2*sqajk);
              t1*=sqajk;
            }
            
            if(flag&0x3 && !zero_force){// calculate forces on ions
              cdouble dEw=pref_eiq*I0*t1*part_jk;
              dEw=(dE-dEw)/ngjki;
              Vector_3 dir=-real(gjki);
              dir.normalize();
              fi[i]+=M1if*real( dEw)*dir;          
            }
            if(flag&(0x8|0x4)){ //electron forces needed
              cdouble dv_ak= (zero_force ?  0.: (djk*gjki)*(v-t1)/(wjk.a*ngjki*ngjki)) +t2;
              cVector_3 dv_bk= zero_force ? cVector_3() :  -gjki*(v-t1)/(2*wjk.a*ngjki*ngjki);
              eterm_deriv(ic1,s1,c1,j1,ic1,s1,c1,k1,M1if*pref_eiq,o,v,dv_ak,dv_ak,dv_bk,dv_bk);
            }
          }
        
          // extra constraint energy
          if(j1==k1 && constraint == HARM && M1e){ 
            cdouble v=conj(wj.a)*wj.a;
            Ew += harm_w0_4 * real(part_jk*v)/wf_norm[s1][c1];
            if(flag&(0x8|0x4)){ //electron forces needed
              cdouble dv_ak=     conj(wj.a);
              cdouble dv_aj_conj =wj.a;
              eterm_deriv(ic1,s1,c1,j1,ic1,s1,c1,k1,harm_w0_4/wf_norm[s1][c1],o,v,dv_aj_conj,dv_ak,cVector_3(),cVector_3());
            }
          }
# if 1
          // second block
          // e-e interaction
          for(int s2=s1;s2<2;s2++){
            int ic2=0; // starting index of the wp for current electron
        
            for(int c2=0 ;c2<ne[s2];ic2+=nspl[s2][c2],c2++){ // incrementing block2 wp address
              if(s1==s2 && c2<=c1)
                continue;

              double pref_ee=coul_pref/(wf_norm[s1][c1]*wf_norm[s2][c2]);
              double dE=0.;
              for(int j2=0;j2<nspl[s2][c2];j2++){
                
                cdouble cj2(split_c[s2][ic2+j2][0],split_c[s2][ic2+j2][1]);
                WavePacket& wj2=wp[s2][ic2+j2];
                
                OverlapDeriv o2;
                if(flag&(0x8|0x4)) //electron forces needed
                  o2.set1(wj2);

                for(int k2=j2;k2<nspl[s2][c2];k2++){
                  int M2a=(j2==k2 ? 1: 2);
                  double M2e, M2f;
                  _mytie(M2e,M2f)=check_part1(s1,ic1+j1,ic1+k1,s2,ic2+j2,ic2+k2);

                  cdouble ck2(split_c[s2][ic2+k2][0],split_c[s2][ic2+k2][1]);

                  WavePacket wk2=wp[s2][ic2+k2];
                  if(pbc)
                    move_to_image(wj2,wk2);
                  WavePacket wjk2=conj(wj2)*wk2;
                  cdouble I02 = wjk2.integral();


                  cdouble part_jk2=conj(cj2)*ck2;
                  
                  cVector_3 djk2=wjk2.b/(2*wjk2.a);
                  cVector_3 ddv=djk-djk2;
                  cdouble dd=ddv.norm();
                  cdouble aa=1./sqrt(1./wjk.a+1./wjk2.a);
                  //double ww1=wj.get_width();
                  //double ww2=wj2.get_width();
                  cdouble v=cerf_div(dd,aa);
                  double pref_eeq=pref_ee*(qe[s1][ic1+j1]+qe[s1][ic1+k1])*(qe[s2][ic2+j2]+qe[s2][ic2+k2])*0.25;
                  cdouble Vj1j2k1k2=pref_eeq*I0*I02*v*part_jk*part_jk2;
                  Eee+=M1e*M2e*real(Vj1j2k1k2);
                  //Eee+=(j1==k1 || j2==k2 ? 1: 2)*real(Vj1j2k1k2);
                  

                  cdouble t1, t2;
                  bool zero_force=(fabs(real(dd))+fabs(imag(dd))<1e-10);
                  if(flag&(0x8|0x4) ){ //electron forces needed
                    cdouble arg=dd*aa;
                    t1=two_over_sqr_pi*exp(-arg*arg)*aa;
                    t2=t1*aa*aa/2;
                  
                    
                    cdouble dv_ak1= (zero_force ?  0.: (djk*ddv)*(v-t1)/(wjk.a*dd*dd)) +t2/(wjk.a*wjk.a);
                    cVector_3 dv_bk1= zero_force ? cVector_3() :  -ddv*(v-t1)/(2*wjk.a*dd*dd);
                    eterm_deriv(ic1,s1,c1,j1,ic1,s1,c1,k1,M1f*M2f*pref_eeq*I02*part_jk2,o,v,dv_ak1,dv_ak1,dv_bk1,dv_bk1);
                    
                    o2.set2(wk2,&I02);
                    cdouble dv_ak2= (zero_force ?  0.: (-djk2*ddv)*(v-t1)/(wjk2.a*dd*dd)) +t2/(wjk2.a*wjk2.a);
                    cVector_3 dv_bk2= zero_force ? cVector_3() :  ddv*(v-t1)/(2*wjk2.a*dd*dd);
                    eterm_deriv(ic2,s2,c2,j2,ic2,s2,c2,k2,M1f*M2f*pref_eeq*I0*part_jk,o2,v,dv_ak2,dv_ak2,dv_bk2,dv_bk2);
                  }
             
                }// k2
              } // j2
              
            } //c2
          }// s2
# endif
        } // k1
      }// j1
      ic1+=nspl[s1][c1]; // incrementing block1 wp address
    }// c1
  } // s1       

  // transforming the forces to physical coordinates
  if(flag&(0x8|0x4) && !(flag&0x10))
    get_el_forces((flag&(~0x4))|0x8,fe_x,fe_p,fe_w,fe_pw,fe_c); // flag change: electronic forces were cleared already

  if(calc_ii)
    interaction_ii(flag,fi);
  
  return 1;
}


///\en Calcualtes the overlap between two electrons taking all split WPs into account.
///    Norms must be pre-calculated.
cdouble  AWPMD_split::overlap(int ic1, int s1, int c1,int ic2, int s2, int c2){
  cdouble sum(0,0);
  for(int j1=0;j1<nspl[s1][c1];j1++){
    double cj_re=split_c[s1][ic1+j1][0];
    double cj_im=split_c[s1][ic1+j1][1];
    cdouble ccj=cdouble(cj_re,-cj_im);
    WavePacket wj=wp[s1][ic1+j1]; 
    for(int k2=0;k2<nspl[s2][c2];k2++){
      double ck_re=split_c[s2][ic2+k2][0];
      double ck_im=split_c[s2][ic2+k2][1];
      cdouble ck=cdouble(ck_re,ck_im);

      WavePacket wk=wp[s2][ic2+k2];
      if(pbc) 
        move_to_image(wj,wk);
      WavePacket wjk=conj(wj)*wk;
      sum+=ccj*ck*wjk.integral();
    }
  }
  return sum/sqrt(wf_norm[s1][c1]*wf_norm[s2][c2]);
}


/// adds the derivatives of Y in the term v*Y[s](c2,c1)
void AWPMD_split::y_deriv(cdouble v,int s,int c2, int c1){
  int ic=0;
  for(int c=0;c<ne[s];c++){
    for(int j=0;j<nspl[s][c];j++){
      E_der[s][ic]+=real((-M[s](c2,ic)*Y[s](c,c1)-Y[s](c2,c)*conj(M[s](c1,ic)))*v);
      ic++;
    }
  }
}

///\en Calculates interaction in the system of ni ions + electrons 
/// the electonic subsystem must be previously setup by set_electrons, ionic by set_ions
/// 0x1   -- give back ion forces \n
/// 0x2   -- add ion forces to the existing set \n
/// 0x4   -- calculate electronic forces  \n
/// 0x8   -- add electronic forces to the existing arrays \n
/// 0x10  -- calculate internal electronic derivatives only: \n
///          will not update electronic force arrays, which may be NULL, \n
///          the forces may be obtained then using \ref get_el_forces() for all WPs \n
///          or separately for each WP using \ref get_wp_force()
/// if PBCs are used the coords must be within a range [0, cell)
int AWPMD_split::interaction(int flag, Vector_3P fi, Vector_3P fe_x, 
                             Vector_3P fe_p, double *fe_w, double *fe_pw, Vector_2P fe_c){
  //if(approx==HARTREE)
    //return interaction_hartree(flag,fi,fe_x,fe_p,fe_w,fe_pw,fe_c);

  //0. resize arrays if needed
  resize(flag);
  //1. clearing forces
  clear_forces(flag,fi,fe_x,fe_p,fe_w,fe_pw,fe_c);
  // calculate block norms and (optionally) derivatives
  calc_norms(flag);  

  //2. calculating overlap matrix
  for(int s=0;s<2;s++){
    int nes = ne[s];
    if(nes == 0) continue;

    int ik=0;
    for(int k=0;k<nes;k++){
      
      Y[s].set(k,k,1.);  // Diagonal elements (=1)
      Oflg[s](k,k) = 1;
      int il=0;
      for(int l=0/*k+1*/;l<nes;il+=nspl[s][l],l++){ // incrementing block2 wp address
        if(l<k+1)
          continue;
        cdouble Okl=overlap(ik,s,k,il,s,l);
        Y[s].set(k,l,Okl);
        Oflg[s](k,l) = norm(Okl) > ovl_tolerance;
      }
      ik+=nspl[s][k]; // incrementing block1 wp address
    }
    O[s] = Y[s];  // save overlap matrix

    //3. inverting the overlap matrix
    int info=0;
    if(nes && approx!=HARTREE){
      /*FILE *f1=fopen(fmt("matrO_%d.d",s),"wt");
      fileout(f1,Y[s],"%15g");
      fclose(f1);8*/
      
      ZPPTRF("L",&nes,Y[s].arr,&info);
      // analyze return code here
      if(info<0)
        return LOGERR(info,fmt("AWPMD.interacton: call to ZPTRF failed (exitcode %d)!",info),LINFO); 
      ZPPTRI("L",&nes,Y[s].arr,&info);
      if(info<0)
        return LOGERR(info,fmt("AWPMD.interacton: call to ZPTRI failed (exitcode %d)!",info),LINFO); 

      
      /*f1=fopen(fmt("matrY_%d.d",s),"wt");
      fileout(f1,Y[s],"%15g");
      fclose(f1);*/
    }

    // Clearing matrices for electronic forces
    /*
    if(flag&0x4){
      Te[s].Set(0.);
      Tei[s].Set(0.);
    }*/
  }
# if 1
  // calculating the M matrix
  if(norm_needed || ( flag&(0x8|0x4) && approx!=HARTREE ) ){
    for(int s1=0;s1<2;s1++){
      if(approx!=HARTREE){
        M[s1].Set(0.);
        L[s1].Set(0.);
        if(norm_needed)
          Norm[s1].Set(0.);
      }
      else
        Lh[s1].assign(nvar[s1],0.);

      int indn=(nvar[s1]/10)*8;
      int ic1=0;
      for(int c1=0;c1<ne[s1];c1++){
        // L[s](c1,c1*)=0
        /*for(int j=0;j<nspl[s1][c1];j++){
          for(int i=0;i<10;i++)
            L[s1](c1,10*(ic1+j)+i)=0.;
        }*/

        int ic2=0;     
        for(int c2=0;c2<ne[s1];ic2+=nspl[s1][c2],c2++){
          if(approx==HARTREE && c1!=c2) // only diagonals for Hartree
            continue;
          if(c2<c1) // taken assymmetry into account
            continue;
          if( !Oflg[s1](c1,c2) ) 
            continue; // non-overlapping WPs

          
          if(approx==HARTREE && norm_needed) // initializing the matrices with zero
            Normh[s1][c1].Set(0);
          
          
          double sq_norm12=sqrt(wf_norm[s1][c1]*wf_norm[s1][c2]);
          double pref=1./(sq_norm12); 
          
         
          // WP blocks:
          for(int j1=0;j1<nspl[s1][c1];j1++){
            cdouble cj1(split_c[s1][ic1+j1][0],split_c[s1][ic1+j1][1]);
            WavePacket wj1=wp[s1][ic1+j1];  

            OverlapDeriv o12;
            o12.set1(wj1);

           for(int k2=(c1==c2 ? j1: 0); k2<nspl[s1][c2];k2++){
              double M12= (c1==c2 && j1==k2) ? 0.5 : 1;

              cdouble ck2(split_c[s1][ic2+k2][0],split_c[s1][ic2+k2][1]);

              // electrons kinetic energy
              WavePacket wk2=wp[s1][ic2+k2];
              if(pbc)
                move_to_image(wj1,wk2);
              
              WavePacket wjk12=conj(wj1)*wk2;
              cdouble I012 = wjk12.integral();
              cdouble part_jk12=conj(cj1)*ck2;     
              
              o12.set2(wk2,&I012);


              cdouble der_k[10], der_j[10];
              // over a_k_re
              der_k[0]=   part_jk12*o12.da2_re();           
              // over a_k_im
              der_k[1]=   part_jk12*o12.da2_im();    
              // over a_j_re
              der_j[0]=   part_jk12*o12.da1_re();           
              // over a_j_im
              der_j[1]=   part_jk12*o12.da1_im();           

              for(int i=0;i<3;i++){
                // over b_k_re
                der_k[2+2*i]=   part_jk12*o12.db2_re(i);           
                // over b_k_im
                der_k[2+2*i+1]= part_jk12*o12.db2_im(i);      
                // over b_j_re
                der_j[2+2*i]=   part_jk12*o12.db1_re(i);                      
                // over b_j_im
                der_j[2+2*i+1]= part_jk12*o12.db1_im(i);             
              }

              // over ck_re
              der_k[8]=conj(cj1)*I012;  
              // over ck_im
              der_k[9]=i_unit*conj(cj1)*I012;
              // over cj_re
              der_j[8]=ck2*I012;  
              // over cj_im
              der_j[9]=-i_unit*ck2*I012;  

              if(j1==0){ // add all together instead of adding by parts
                cdouble t=-O[s1](c1,c2)/pref;   //-part_jk12*I012; 
                for(int i=0;i<8;i++){
                  der_j[i]+=t*wf_norm_der[s1][8*(ic1+j1)+i];
                  der_k[i]+=t*wf_norm_der[s1][8*(ic2+k2)+i];
                }
                der_j[8]+=t*wf_norm_der[s1][indn+2*(c1+j1)];
                der_j[9]+=t*wf_norm_der[s1][indn+2*(c1+j1)+1];
                der_k[8]+=t*wf_norm_der[s1][indn+2*(c2+k2)];
                der_k[9]+=t*wf_norm_der[s1][indn+2*(c2+k2)+1];
              }
            
              if(approx!=HARTREE){
                for(int i=0;i<10;i++){
                  L[s1](c1,10*(ic2+k2)+i)+=M12*pref*der_k[i];
                  L[s1](c2,10*(ic1+j1)+i)+=M12*pref*conj(der_j[i]);
                }

                // M=Y*L
                for(int g=0;g<ne[s1];g++){
                  for(int i=0;i<10;i++){
                    M[s1](g,10*(ic2+k2)+i)+=M12*pref*Y[s1](g,c1)*der_k[i];
                    M[s1](g,10*(ic1+j1)+i)+=M12*pref*Y[s1](g,c2)*conj(der_j[i]);
                  }
                }//g
              }
              else{ // HARTREE
                for(int i=0;i<10;i++){
                  Lh[s1][10*(ic2+k2)+i]+=M12*pref*der_k[i];
                  Lh[s1][10*(ic1+j1)+i]+=M12*pref*der_k[i];
                }
              }
              if(norm_needed){ // filling part of norm matrix
                o12.calc_der_overlap(false,conj(cj1),ck2);
                if(approx!=HARTREE){
                  for(int i=0;i<10;i++)  // 10x10
                    for(int j=0;j<10;j++)
                      Norm[s1](10*(ic1+j1)+i,10*(ic2+k2)+j)=-2*imag(o12.IDD((int)i,(int)j))/one_h;              
                }
                else{
                  for(int i=0;i<10;i++)  // 10x10
                    for(int j=0;j<10;j++)
                      Normh[s1][c1](10*j1+i,10*k2+j)=-2*imag(o12.IDD((int)i,(int)j))/one_h;              
                }
              }
            } // k2
          } //j1
        }// c2
        ic1+=nspl[s1][c1]; // incrementing block1 wp address
      }// c1
      if(norm_needed){  // filling the rest of norm_matrix
        if(approx!=HARTREE){
          int ic1=0;
          for(int c1=0;c1<ne[s1];c1++){
            int ic2=0;     
            for(int c2=0;c2<ne[s1];ic2+=nspl[s1][c2],c2++){
              if(c2<c1)
                continue;
              for(int j1=0;j1<10*nspl[s1][c1];j1++){
                for(int k2=(c1==c2 ? j1: 0); k2<10*nspl[s1][c2];k2++){
                  cdouble y=0.;
                  for(int g=0;g<ne[s1];g++)
                    y+=conj(L[s1](g,10*ic1+j1))*M[s1](g,10*ic2+k2);
                  cdouble elm=-conj(L[s1](c2,10*ic1+j1))*wf_norm_der[s1][10*ic2+k2]-L[s1](c1,10*ic2+k2)*wf_norm_der[s1][10*ic1+j1]-O[s1](c1,c2)*wf_norm_der[s1][10*ic2+k2]*wf_norm_der[s1][10*ic1+j1];
                  elm-=y;
                  elm*=Y[s1](c2,c1);
                  Norm[s1](10*ic1+j1,10*ic2+k2)+=-2*imag(elm)/one_h;
                }// k2
              } // j1
            }// c2
            ic1+=nspl[s1][c1]; // incrementing block1 wp address
          }
        }
        else{ // HARTREE
          int ic1=0;
          for(int c1=0;c1<ne[s1];c1++){
            for(int j1=0;j1<10*nspl[s1][c1];j1++){
              for(int k2=j1; k2<10*nspl[s1][c1];k2++){
                cdouble y=conj(Lh[s1][10*ic1+j1])*Lh[s1][10*ic1+k2];
                cdouble elm=-conj(Lh[s1][10*ic1+j1])*wf_norm_der[s1][10*ic1+k2]-Lh[s1][10*ic1+k2]*wf_norm_der[s1][10*ic1+j1]-wf_norm_der[s1][10*ic1+k2]*wf_norm_der[s1][10*ic1+j1];
                elm-=y;
                Normh[s1][c1](j1,k2)+=-2*imag(elm)/one_h;
              }
            }
            ic1+=nspl[s1][c1]; // incrementing block1 wp address
          }
        }
      } // norm 
    } // s1    
  } // flag
# endif
  Vector_3 ndr;
  //  calculating single particle contribution
  Edc = Edk = Eee = Ew = 0.;
  for(int s1=0;s1<2;s1++){
    Ee[s1]=Eei[s1]=0.;
    int ic1=0;
    for(int c1=0;c1<ne[s1];c1++){
      
      //double Ee1=0., Ew1=0., Eei1=0.;
      
      int ic2=0;
      for(int c2=0;c2<ne[s1];ic2+=nspl[s1][c2],c2++){
        if( !Oflg[s1](c1,c2) ) 
            continue; // non-overlapping WPs

        if(c2<c1) // taken as M factor into account
          continue;
       
        double sq_norm12=sqrt(wf_norm[s1][c1]*wf_norm[s1][c2]);
        double pref=-h2_me/(2*sq_norm12); // ekin
        double pref_ei=coul_pref/sq_norm12;

        cdouble yy;
        if(approx==HARTREE)
          yy=1.;
        else 
          yy=Y[s1](c2,c1);
        cdouble Tc1c2=0.;

        if((approx!=HARTREE || c2==c1) && Oflg[s1](c1,c2)){ // only diagonal terms for Hartree, and overlap nonzero
          // WP blocks:
          for(int j1=0;j1<nspl[s1][c1];j1++){
            cdouble cj1(split_c[s1][ic1+j1][0],split_c[s1][ic1+j1][1]);
            WavePacket wj1=wp[s1][ic1+j1];  

            OverlapDeriv o12;
            if(flag&(0x8|0x4)) //electron forces needed
              o12.set1(wj1);

            for(int k2=(c1==c2 ? j1: 0); k2<nspl[s1][c2];k2++){
              int M12=(c1==c2 && j1==k2 ? 1: 2);
              double M12pe, M12pf;
              _mytie(M12pe,M12pf)=check_part1(s1,ic1+j1,ic2+k2)*M12;
              
                  

              cdouble ck2(split_c[s1][ic2+k2][0],split_c[s1][ic2+k2][1]);
            

              // electrons kinetic energy
              WavePacket wk2=wp[s1][ic2+k2];
              if(pbc)
                move_to_image(wj1,wk2);
              
              WavePacket wjk12=conj(wj1)*wk2;
              cdouble I012 = wjk12.integral();
              cdouble part_jk12=conj(cj1)*ck2;           
              cVector_3 djk12=wjk12.b/(2.*wjk12.a);

              // kinetic energy contribution
              if(M12pf){
                cVector_3 v1=conj(wj1.b)*wk2.a-wk2.b*conj(wj1.a);
                cdouble v=(v1*v1)/wjk12.a;
                v-=6.*wk2.a*conj(wj1.a);
                v/=wjk12.a;
                //v=1.; 
         
                double dE=M12pe*real(part_jk12*I012*v*yy)*pref;  
                //double dE=real(yy);
                //Tc1c2+=1.;
                Ee[s1] += dE;
                Eep[s1][ic1+j1]+= 0.5*dE; //per particle energy
                Eep[s1][ic2+k2]+= 0.5*dE;

                // diagonal term
                if(M12==1)
                  Edk+=dE;
         
                if(flag&(0x8|0x4)){ //electron forces needed
                  cVector_3 tv=wk2.b*conj(wj1.a)-conj(wj1.b)*wk2.a;
                  cdouble ajk2=wjk12.a*wjk12.a;
                  cdouble ajk3=ajk2*wjk12.a;
                  cdouble dv_aj_conj=     -2*wk2.a*(3*wjk12.a*wk2.a-tv*wjk12.b)/ajk3;
                  cdouble dv_ak=-2*conj(wj1.a)*(3*wjk12.a*conj(wj1.a)+tv*wjk12.b)/ajk3;
                  cVector_3 dv_bj_conj=-2*tv*wk2.a/ajk2;
                  cVector_3 dv_bk=2*tv*conj(wj1.a)/ajk2;

                  
                  o12.set2(wk2,&I012);
                  // calculate full derivative of the term pref*conj(cj)*ck*Ijk*vjk/sqrt(nrm(s1)*nrm(s2))
                  // denominator must be included in pref
                  eterm_deriv(ic1,s1,c1,j1,ic2,s1,c2,k2,M12pf*pref*yy,o12,v,dv_aj_conj,dv_ak,dv_bj_conj,dv_bk);

                  Tc1c2+=M12pf*part_jk12*I012*v*pref; // all without Y
                }
              } // M12pe
              // e-i energy
              cdouble sum=0.;
              for(int i=0;i<ni;i++){  // ions loop 
                double M12pie, M12pif;
                _mytie(M12pie,M12pif)=check_part1ei(s1,ic1+j1,ic2+k2,i)*M12;
                if(!M12pif)
                  continue;
                cVector_3 gjki=djk12 - cVector_3(xi[i]);
                
                if(pbc) // correcting the real part (distance) according to PBC
                  gjki=rcell1(gjki,cell,pbc);
                  //-Igor- gkli=cVector_3(real(gkli).rcell1(cell,pbc),imag(gkli));
                
                cdouble ngjki=gjki.norm();
                cdouble sqajk=sqrt(wjk12.a);
                //cdouble ttt = cerf_div(ngjki,c);
                cdouble v=cerf_div(ngjki,sqajk);
                double pref_eiq=pref_ei*qi[i]*(qe[s1][ic1+j1]+qe[s1][ic2+k2])*0.5;
                cdouble dE=pref_eiq*I012*v*part_jk12;
                sum+=M12pie*dE; // sum without Y
                dE*=yy;

                Eeip[s1][ic1+j1]+= 0.25*real(dE); //per particle energy
                Eeip[s1][ic2+k2]+= 0.25*real(dE);
                Eiep[i]+=0.5*real(dE);

                cdouble t1, t2;
                bool zero_force=(fabs(real(ngjki))+fabs(imag(ngjki))<1e-10);
                if(flag&(0x8|0x4|0x3) ){ //some derivatives needed
                  cdouble arg=ngjki*sqajk;
                  t1=two_over_sqr_pi*exp(-arg*arg);
                  t2=t1/(2*sqajk);
                  t1*=sqajk;
                }
                
                if(flag&0x3 && !zero_force){// calculate forces on ions
                  cdouble dEw=(pref_eiq*yy)*I012*t1*part_jk12;
                  dEw=(dE-dEw)/ngjki;
                  Vector_3 dir=-real(gjki);
                  dir.normalize();
                  fi[i]+=M12pif*real( dEw)*dir;          
                }
                if(flag&(0x8|0x4)){ //electron forces needed
                  cdouble dv_ak= (zero_force ?  0.: (djk12*gjki)*(v-t1)/(wjk12.a*ngjki*ngjki)) +t2;
                  cVector_3 dv_bk= zero_force ? cVector_3() :  -gjki*(v-t1)/(2*wjk12.a*ngjki*ngjki);
                  eterm_deriv(ic1,s1,c1,j1,ic2,s1,c2,k2,M12pif*pref_eiq*yy,o12,v,dv_ak,dv_ak,dv_bk,dv_bk);
                }
              }
              // ions energy contribution
              Eei[s1] += real(sum*yy);

              // diagonal term
              if(M12==1)
                Edc+= real(sum*yy);;

              if(flag&(0x8|0x4)) //electron forces needed
                Tc1c2+=sum;

              // extra constraint energy, will only arrive here when M12p=1
              if((c1==c2 && j1==k2) && constraint == HARM){ 
                cdouble v=conj(wj1.a)*wk2.a;
                double dE= M12pe * harm_w0_4 * real(part_jk12*v)/wf_norm[s1][c1];
                Ew += dE;

                Ewp[s1][ic1+j1]+= dE; //per particle energy
                

                if(flag&(0x8|0x4)){ //electron forces needed
                  cdouble dv_ak=     conj(wj1.a);
                  cdouble dv_aj_conj =wk2.a;
                  eterm_deriv(ic1,s1,c1,j1,ic1,s1,c1,k2,M12pf*harm_w0_4/wf_norm[s1][c1],o12,v,dv_aj_conj,dv_ak,cVector_3(),cVector_3());
                }
              }
            }// k2
          }// j1
          if(flag&(0x8|0x4) && approx!=HARTREE){ //electron forces needed
            // adding Y derivative term for all variables (te and tei terms)
            y_deriv(Tc1c2,s1,c2,c1);
          }
        } // !HARTREE or spins or overlap
      


# if 1 // pair by pair sum
        // second block
        // e-e interaction
        for(int s2=s1;s2<2;s2++){
          if(approx==HARTREE && c1!=c2) // only Vkmkm terms for Hartree
               continue;
          if(/*s1==s2 &&*/ c2<c1) // pair selection term
             continue;
                  
          int ic3=0; // starting index of the wp for current electron
          for(int c3=0 ;c3<ne[s2];ic3+=nspl[s2][c3],c3++){ // incrementing block2 wp address
            if(approx==HARTREE && s1==s2 && c3==c1) // [Vkkmn contribution for same spin is 0], also  no Vkkkk terms for Hartree (=)
              continue;
            if(s1==s2 && c3<=c1) // and general Coulomb sum: i<j (<) //??
              continue;
            
            int ic4=0;   
            for(int c4=0; c4<ne[s2];ic4+=nspl[s2][c4],c4++){
              if(approx==HARTREE && c4!=c3) // only Vkmkm terms for Hartree
                continue;
              //if(s1==s2 && c4==c2) // Vklmm contribution for same spin is 0 for antisymmetrized approximations, also Vmmmm is 0 for Hartree
              //  continue;
              if(s1==s2 && c4<=c2) // pair selection term: Vklmm
                continue;
              if(/*s1==s2 &&*/ c4<c3) // pair selection term
                continue;

              double sq_norm34=sqrt(wf_norm[s2][c3]*wf_norm[s2][c4]);
              double pref_ee=coul_pref/(sq_norm12*sq_norm34);
                  
                
              cdouble yy;
              double K=1.;
              if(approx==HARTREE){
                yy=1.;
              }
              else{
                if(s1==s2){ // same spin antisymmetrized term
                  yy=Y[s1](c2,c1)*Y[s1](c4,c3)-Y[s1](c2,c3)*Y[s1](c4,c1);  // check the order of c
                }
                else{ // different spin atisymmetrized term
                  //if(approx==UHF)
                    //K=2.;
                  yy=K*Y[s2](c4,c3)*Y[s1](c2,c1);                 
                }
                //yy=1.;
                //yy=Y[s1](c2,c1)*Y[s1](c4,c3);
                //yy=Y[s1](c2,c3)*Y[s1](c4,c1);
                //yy=Y[s1](c2,c1)*Y[s1](c4,c3)-Y[s1](c2,c3)*Y[s1](c4,c1);
              }
          
              cdouble Vy=0.;
              // WP blocks: Vc1c3c2c4
              for(int j1=0;j1<nspl[s1][c1];j1++){
                cdouble cj1(split_c[s1][ic1+j1][0],split_c[s1][ic1+j1][1]);
                WavePacket wj1=wp[s1][ic1+j1];  

                OverlapDeriv o12,o14;
                if(flag&(0x8|0x4)){ //electron forces needed
                  o12.set1(wj1);
                  o14.set1(wj1);
                }

                for(int k2=(approx==HARTREE ? j1: 0); k2<nspl[s1][c2];k2++){
                  int M12=(c1==c2 && j1==k2 ? 1: 2);
                  
                 
                  cdouble ck2(split_c[s1][ic2+k2][0],split_c[s1][ic2+k2][1]);
                

                  WavePacket wk2=wp[s1][ic2+k2];
                  if(pbc)
                    move_to_image(wj1,wk2);
                  
                  WavePacket wjk12=conj(wj1)*wk2;
                  cdouble I012 = wjk12.integral();
                  cdouble part_jk12=conj(cj1)*ck2;           
                  cVector_3 djk12=wjk12.b/(2.*wjk12.a);

                  if(flag&(0x8|0x4)) //electron forces needed
                    o12.set2(wk2,&I012);

                  
                  for(int j3=0;j3<nspl[s2][c3];j3++){
                    
                    cdouble cj3(split_c[s2][ic3+j3][0],split_c[s2][ic3+j3][1]);
                    WavePacket& wj3=wp[s2][ic3+j3];
                    
                    OverlapDeriv o34, o32;
                    if(flag&(0x8|0x4)){ //electron forces needed
                      o34.set1(wj3);
                      o32.set1(wj3);
                      //o32.set2(wk2);
                    }
                      
                    // 3-2
                    WavePacket wjk32;
                    cdouble I032=0., part_jk32;
                    cVector_3 djk32;
                    if(s1==s2 || approx==UHF){
                      WavePacket wk2=wp[s2][ic2+k2];
                      if(pbc)
                         move_to_image(wj3,wk2);
                      wjk32=conj(wj3)*wk2;
                      I032 = wjk32.integral();
                      part_jk32=conj(cj3)*ck2;
                      djk32=wjk32.b/(2*wjk32.a);
                      if(flag&(0x8|0x4)) //electron forces needed
                        o32.set2(wk2,&I032);
                    }
                    

                    for(int k4=(approx==HARTREE ? j3: 0);k4<nspl[s2][c4];k4++){
                      int M34=(c3==c4 && j3==k4 ? 1: 2);
                       
                      

                      double M0;
                      if(approx==HARTREE)
                        M0=M12*M34;
                      else{
                        M0= (c1==c2 && c3==c4 ? 1: 2); // will have exchange term for different pairs instead of M12*M34 factor
                      }
                      double Me, Mf;
                      _mytie(Me,Mf)=check_part1(s1,ic1+j1,ic2+k2,s2,ic3+j3,ic4+k4)*M0;
                      if(!Mf)
                        continue;

                      cdouble ck4(split_c[s2][ic4+k4][0],split_c[s2][ic4+k4][1]);

                      // 3-4
                      WavePacket wk4=wp[s2][ic4+k4];
# if 1
                      if(pbc)
                        move_to_image(wj3,wk4);
                      WavePacket wjk34=conj(wj3)*wk4;
                      cdouble I034 = wjk34.integral();
                      cdouble part_jk34=conj(cj3)*ck4;
                      cVector_3 djk34=wjk34.b/(2*wjk34.a);
                      
                      cVector_3 ddv=djk12-djk34;
                      cdouble dd=ddv.norm();
                      cdouble aa=1./sqrt(1./wjk12.a+1./wjk34.a);
                      cdouble v=cerf_div(dd,aa);
                      double pref_eeq=pref_ee*(qe[s1][ic1+j1]+qe[s1][ic2+k2])*(qe[s2][ic3+j3]+qe[s2][ic4+k4])*0.25;
                      cdouble Vj1j3k2k4=pref_eeq*I012*I034*v*part_jk12*part_jk34;  // Vklmn
                      double dE=Me*real(Vj1j3k2k4*yy);
                      Eee+=dE;
                      
                      Eeep[s1][ic1+j1]+=0.25*dE; // per-particle energy
                      Eeep[s1][ic2+k2]+=0.25*dE;
                      Eeep[s2][ic3+j3]+=0.25*dE;
                      Eeep[s2][ic4+k4]+=0.25*dE;

                      // diagonal term
                      if(M12==1 && M34 ==1)
                        Edc+=dE;

                      bool zero_force=(fabs(real(dd))+fabs(imag(dd))<1e-10);
                      if(flag&(0x8|0x4) ){ //electron forces needed
                        cdouble arg=dd*aa;
                        cdouble t1=two_over_sqr_pi*exp(-arg*arg)*aa;
                        cdouble t2=t1*aa*aa/2;
                      
                        
                        cdouble dv_ak1= (zero_force ?  0.: (djk12*ddv)*(v-t1)/(wjk12.a*dd*dd)) +t2/(wjk12.a*wjk12.a);
                        cVector_3 dv_bk1= zero_force ? cVector_3() :  -ddv*(v-t1)/(2*wjk12.a*dd*dd);
                        eterm_deriv(ic1,s1,c1,j1,ic2,s1,c2,k2,Mf*pref_eeq*I034*part_jk34*yy,o12,v,dv_ak1,dv_ak1,dv_bk1,dv_bk1);
                        
                        o34.set2(wk4,&I034);
                        cdouble dv_ak2= (zero_force ?  0.: (-djk34*ddv)*(v-t1)/(wjk34.a*dd*dd)) +t2/(wjk34.a*wjk34.a);
                        cVector_3 dv_bk2= zero_force ? cVector_3() :  ddv*(v-t1)/(2*wjk34.a*dd*dd);
                        eterm_deriv(ic3,s2,c3,j3,ic4,s2,c4,k4,Mf*pref_eeq*I012*part_jk12*yy,o34,v,dv_ak2,dv_ak2,dv_bk2,dv_bk2);

                        Vy+=Me*Vj1j3k2k4; // all without yy
                      }// force
# endif
# if 1
                      // 1-4
                      if(approx!=HARTREE && (approx==UHF || s1==s2)){
                        wk4=wp[s2][ic4+k4];
                        if(pbc)
                          move_to_image(wj1,wk4);
                        WavePacket wjk14=conj(wj1)*wk4;
                        cdouble I014 = wjk14.integral();
                        cdouble part_jk14=conj(cj1)*ck4;
                        cVector_3 djk14=wjk14.b/(2*wjk14.a);
                        
                        cVector_3 ddv=djk32-djk14;
                        cdouble dd=ddv.norm();
                        cdouble aa=1./sqrt(1./wjk32.a+1./wjk14.a);
                        cdouble v=cerf_div(dd,aa);

                        
                        double pref_eeq=pref_ee*(qe[s1][ic1+j1]+qe[s2][ic4+k4])*(qe[s2][ic3+j3]+qe[s1][ic2+k2])*0.25;
                        cdouble Vj1j3k4k2=pref_eeq*I032*I014*v*part_jk32*part_jk14;  // Vklmn
                        double dE=-Me*real(Vj1j3k4k2*yy);
                        Eee+=dE;
                        
                        Eeep[s1][ic1+j1]+=0.25*dE; // per-particle energy
                        Eeep[s1][ic2+k2]+=0.25*dE;
                        Eeep[s2][ic3+j3]+=0.25*dE;
                        Eeep[s2][ic4+k4]+=0.25*dE;


                        bool zero_force=(fabs(real(dd))+fabs(imag(dd))<1e-10);
                        if(flag&(0x8|0x4) ){ //electron forces needed
                          cdouble arg=dd*aa;
                          cdouble t1=two_over_sqr_pi*exp(-arg*arg)*aa;
                          cdouble t2=t1*aa*aa/2;
                        
                          o14.set2(wk4,&I014);
                          cdouble dv_ak1= (zero_force ?  0.: (-djk14*ddv)*(v-t1)/(wjk14.a*dd*dd)) +t2/(wjk14.a*wjk14.a);
                          cVector_3 dv_bk1= zero_force ? cVector_3() :  ddv*(v-t1)/(2*wjk14.a*dd*dd);
                          eterm_deriv(ic1,s1,c1,j1,ic4,s2,c4,k4,-Mf*pref_eeq*I032*part_jk32*yy,o14,v,dv_ak1,dv_ak1,dv_bk1,dv_bk1);
                          
                          
                          cdouble dv_ak2= (zero_force ?  0.: (djk32*ddv)*(v-t1)/(wjk32.a*dd*dd)) +t2/(wjk32.a*wjk32.a);
                          cVector_3 dv_bk2= zero_force ? cVector_3() :  -ddv*(v-t1)/(2*wjk32.a*dd*dd);
                          eterm_deriv(ic3,s2,c3,j3,ic2,s1,c2,k2,-Mf*pref_eeq*I014*part_jk14*yy,o32,v,dv_ak2,dv_ak2,dv_bk2,dv_bk2);

                          Vy+=-Me*Vj1j3k4k2; // all without yy
                        }// force
                      } // s1==s2
# endif
                    }// k4
                  }// j3
                } // k2  
              } // j1

              // derivative of yy term
              if(approx!=HARTREE && flag&(0x8|0x4)){ //electron forces needed
                //yy=Y[s1](c2,c1)*Y[s1](c4,c3)-Y[s1](c2,c3)*Y[s1](c4,c1); 
# if 1
                if(s1==s2){
                  // yy
                  y_deriv(Vy*Y[s1](c4,c3),s1,c2,c1);
                  y_deriv(Vy*Y[s1](c2,c1),s1,c4,c3);
                  y_deriv(-Vy*Y[s1](c4,c1),s1,c2,c3);
                  y_deriv(-Vy*Y[s1](c2,c3),s1,c4,c1);                  
                }
                else{
                  y_deriv(K*Vy*Y[s2](c4,c3),s1,c2,c1);
                  y_deriv(K*Vy*Y[s1](c2,c1),s2,c4,c3);
                }
# endif
               // y_deriv(Vy*Y[s1](c4,c1),s1,c2,c3);
               // y_deriv(Vy*Y[s1](c2,c3),s1,c4,c1);
              }
            } // c4
          } // c3 
        } // s2
# endif
      }// c2
      ic1+=nspl[s1][c1]; // incrementing block1 wp address
    }// c1
  } // s1       
  
  // transforming the forces to physical coordinates
  if(flag&(0x8|0x4) && !(flag&0x10))
    get_el_forces((flag&(~0x4))|0x8,fe_x,fe_p,fe_w,fe_pw,fe_c); // flag change: electronic forces were cleared already
  
  if(calc_ii)
    interaction_ii(flag,fi);

  return 1;
}