
module mod_force
  use global
  use omp_lib
  use,intrinsic :: ieee_arithmetic

contains



 subroutine lj_cut_coul_dsf_nonewton_SP(step,printflag)
    implicit none
    real*4 :: force,forcelj,forcecoul
    real*4 :: x1,y1,z1,x2,y2,z2
    real*4 :: dx,dy,dz,dr,dr2,dr2i,dr6i,dr12i
    double precision ::  ffx,ffy,ffz
    real*4 :: fudge_factor
    real*4 :: qtmp,dri,r,prefactor,erfcc,erfcd,t
    real*4 :: boxdx,boxdy,boxdz
    real*4 :: scalelj, scalecoul
    integer :: i,j,l,step
    integer :: ioffset
    integer :: itype,jtype,neigh
    integer :: tid,num
    integer :: offset
    integer :: printflag
    integer :: T1,T2,clock_rate,clock_max

    call system_clock(T1,clock_rate,clock_max)


    !dir$ offload begin target(mic:0) in(position: alloc_if(.false.),free_if(.false.)),&
    !dir$ inout(ff: alloc_if(.false.), free_if(.false.)),&
    !dir$ in(nlist: alloc_if(.false.) free_if(.false.)),&
    !dir$ in(nnfirst: alloc_if(.false.) free_if(.false.)),&
    !dir$ in(nnlast: alloc_if(.false.) free_if(.false.)),&
    !dir$ in(q: alloc_if(.false.) free_if(.false.)),&
    !dir$ nocopy(lj1: alloc_if(.false.) free_if(.false.)),&
    !dir$ nocopy(lj2: alloc_if(.false.) free_if(.false.)),&
    !dir$ nocopy(lj3: alloc_if(.false.) free_if(.false.)),&
    !dir$ nocopy(lj4: alloc_if(.false.) free_if(.false.))

    !$omp parallel do schedule(dynamic) reduction(+:potential,e_coul,ffx,ffy,ffz) default(private),&
    !$omp& shared(box,ibox,rcut2,cut_coulsq,alpha,EWALD_P,MY_PIS,np,numAtomType),&
    !$omp& shared(A1,A2,A3,A4,A5),&
    !$omp& shared(e_shift,f_shift),&
    !$omp& shared(position,ff,nlist,nnfirst,nnlast,q,lj1,lj2,lj3,lj4)
    do i = 1 ,np
       x1 = position(i)%x; y1 = position(i)%y; z1 = position(i)%z; itype = position(i)%type
       qtmp = q(i)

       ioffset = numAtomType*(itype-1)
       ffx = 0.0d0; ffy = 0.0d0; ffz = 0.0d0

       !dir$ simd reduction(+:potential,e_coul,ffx,ffy,ffz)
       do j= nnfirst(i),nnlast(i)
          
          neigh = nlist(j)
          jtype = position(neigh)%type

          dx = x1-position(neigh)%x; dy = y1-position(neigh)%y; dz = z1-position(neigh)%z
         
          boxdx = dx*ibox; boxdy = dy*ibox; boxdz = dz*ibox
          boxdx = (boxdx+sign(1/(epsilon(boxdx)),boxdx)) -sign(1/epsilon(boxdx),dx)
          boxdy = (boxdy+sign(1/(epsilon(boxdy)),boxdy)) -sign(1/epsilon(boxdy),dy)
          boxdz = (boxdz+sign(1/(epsilon(boxdz)),boxdz)) -sign(1/epsilon(boxdz),dz)

          dx = dx-box*boxdx; dy = dy-box*boxdy; dz = dz-box*boxdz

          dr2 = dx*dx + dy*dy + dz*dz

          !----lennard jones calculation
          dr2i = 1.0d0/dr2
          dr6i = dr2i*dr2i*dr2i
          offset  = ioffset+jtype
          
          if(dr2.gt.rcut2)dr6i=0.0d0
          forcelj = dr6i*(lj1(offset)*dr6i-lj2(offset))
          potential = potential + dr6i*(dr6i*lj3(offset)-lj4(offset))      


          !---electrostatic calculations
          r = sqrt(dr2)
          dri = 1.0/r
          prefactor =  qtmp*q(neigh)*dri
          erfcd = exp(-alpha*alpha*dr2)
          t = 1.0/(1.0 + EWALD_P*alpha*r)
          erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd

          if(dr2.gt.cut_coulsq)prefactor = 0.0
          forcecoul = prefactor * (erfcc*dri + 2.0*alpha*MY_PIS_INV * erfcd +&
               r*f_shift) * r
          
          e_coul = e_coul + prefactor*(erfcc-r*e_shift-dr2*f_shift)                     
          force   = (forcecoul+forcelj)*dr2i


          ffx = ffx + dx*force
          ffy = ffy + dy*force
          ffz = ffz + dz*force
          
       enddo
       ff(i)%x = ffx; ff(i)%y = ffy; ff(i)%z = ffz
    enddo
    !$omp end parallel do 
    !dir$ end offload

    potential = 0.50d0*potential
    e_coul = 0.50d0*e_coul

    call system_clock(T2,clock_rate,clock_max)

    if(printflag.eq.1)then
       print*,'elapsed time w/o newton SP',real(T2-T1)/real(clock_rate)
       print*,'calculated potentials:',potential*reps, e_coul*reps
    endif
    
  end subroutine lj_cut_coul_dsf_nonewton_SP





 subroutine lj_cut_coul_dsf_nonewton_DP(step,printflag)
    implicit none
    double precision :: force,forcelj,forcecoul
    double precision :: x1,y1,z1,x2,y2,z2
    double precision :: dx,dy,dz,dr,dr2,dr2i,dr6i,dr12i,dri
    double precision :: ffx,ffy,ffz
    double precision :: fudge_factor
    double precision :: qtmp,r,prefactor,erfcc,erfcd,t
    double precision :: boxdx,boxdy,boxdz
    double precision :: scalelj, scalecoul
    integer :: i,j,l,step,ioffset
    integer :: itype,jtype,neigh
    integer :: tid,num
    integer :: offset
    integer :: printflag
    integer :: T1,T2,clock_rate,clock_max

    call system_clock(T1,clock_rate,clock_max)

    !dir$ offload begin target(mic:0) in(pdp: alloc_if(.false.),free_if(.false.)),&
    !dir$ inout(ffdp: alloc_if(.false.), free_if(.false.)),&
    !dir$ in(nlist: alloc_if(.false.) free_if(.false.)),&
    !dir$ in(nnfirst: alloc_if(.false.) free_if(.false.)),&
    !dir$ in(nnlast: alloc_if(.false.) free_if(.false.)),&
    !dir$ in(qdp: alloc_if(.false.) free_if(.false.)),&
    !dir$ nocopy(lj1dp: alloc_if(.false.) free_if(.false.)),&
    !dir$ nocopy(lj2dp: alloc_if(.false.) free_if(.false.)),&
    !dir$ nocopy(lj3dp: alloc_if(.false.) free_if(.false.)),&
    !dir$ nocopy(lj4dp: alloc_if(.false.) free_if(.false.))

    !$omp parallel do schedule(dynamic) reduction(+:potential,e_coul,ffx,ffy,ffz) default(private),&
    !$omp& shared(boxdp,iboxdp,rcut2,cut_coulsq,alphadp,EWALD_Pdp)&
    !$omp& shared(MY_PISdp,MY_PIS_INVdp,np,numAtomType),&
    !$omp& shared(A1dp,A2dp,A3dp,A4dp,A5dp),&
    !$omp& shared(e_shiftdp,f_shiftdp),&
    !$omp& shared(pdp,ffdp,nlist,nnfirst,nnlast,qdp,lj1dp,lj2dp,lj3dp,lj4dp)
    do i = 1 ,np
       x1 = pdp(i)%x; y1 = pdp(i)%y; z1 = pdp(i)%z; itype = pdp(i)%type
       qtmp = qdp(i)

       ioffset = numAtomType*(itype-1)
       ffx = 0.0d0; ffy = 0.0d0; ffz = 0.0d0

       !dir$ simd reduction(+:potential,e_coul,ffx,ffy,ffz)
       do j= nnfirst(i),nnlast(i)
          
          neigh = nlist(j)
          jtype = pdp(neigh)%type

          dx = x1-pdp(neigh)%x; dy = y1-pdp(neigh)%y; dz = z1-pdp(neigh)%z
         
          boxdx = dx*iboxdp; boxdy = dy*iboxdp; boxdz = dz*iboxdp
          boxdx = (boxdx+sign(1/(epsilon(boxdx)),boxdx)) -sign(1/epsilon(boxdx),dx)
          boxdy = (boxdy+sign(1/(epsilon(boxdy)),boxdy)) -sign(1/epsilon(boxdy),dy)
          boxdz = (boxdz+sign(1/(epsilon(boxdz)),boxdz)) -sign(1/epsilon(boxdz),dz)

          dx = dx-boxdp*boxdx; dy = dy-boxdp*boxdy; dz = dz-boxdp*boxdz

          dr2 = dx*dx + dy*dy + dz*dz

          !----lennard jones calculation
          dr2i = 1.0d0/dr2
          dr6i = dr2i*dr2i*dr2i
          offset  = ioffset+jtype
          
          if(dr2.gt.rcut2)dr6i=0.0d0
          forcelj = dr6i*(lj1dp(offset)*dr6i-lj2dp(offset))
          potential = potential + dr6i*(dr6i*lj3dp(offset)-lj4dp(offset))      


          !---electrostatic calculations
          r = sqrt(dr2)
          dri = 1.0/r
          prefactor =  qtmp*qdp(neigh)*dri
          erfcd = exp(-alphadp*alphadp*dr2)
          t = 1.0/(1.0 + EWALD_Pdp*alphadp*r)
          erfcc = t * (A1dp+t*(A2dp+t*(A3dp+t*(A4dp+t*A5dp)))) * erfcd

          if(dr2.gt.cut_coulsq)prefactor = 0.0
          forcecoul = prefactor * (erfcc*dri + 2.0*alpha*MY_PIS_INVdp * erfcd +&
               r*f_shift) * r
          
          e_coul = e_coul + prefactor*(erfcc-r*e_shiftdp-dr2*f_shiftdp)                     
          force   = (forcecoul+forcelj)*dr2i


          ffx = ffx + dx*force
          ffy = ffy + dy*force
          ffz = ffz + dz*force
          
       enddo
       ffdp(i)%x = ffx; ffdp(i)%y = ffy; ffdp(i)%z = ffz
    enddo
    !$omp end parallel do 
    !dir$ end offload
    potential = 0.50d0*potential
    e_coul = 0.50d0*e_coul

    call system_clock(T2,clock_rate,clock_max)

    if(printflag.eq.1)then
       print*,'elapsed time w/o newton DP',real(T2-T1)/real(clock_rate)
       print*,'calculated potentials:',potential*reps, e_coul*reps
    endif
    
  end subroutine lj_cut_coul_dsf_nonewton_DP







  subroutine zero_values(step)
    implicit none
    integer :: step
    integer :: i

    do i =1 ,np
       ff(i)%x = 0.0d0
       ff(i)%y = 0.0d0
       ff(i)%z = 0.0d0

       ffdp(i)%x = 0.0d0
       ffdp(i)%y = 0.0d0
       ffdp(i)%z =0.0d0
    end do

    potential = 0.0d0
    e_coul = 0.0d0



  end subroutine zero_values


  subroutine warm_up_mic_SP
    implicit none




    !dir$ offload_transfer target(mic:0) in(position: alloc_if(.true.) free_if(.false.)),&
    !dir$ in(ff: alloc_if(.true.) free_if(.false.)),&
    !dir$ in(nlist: alloc_if(.true.) free_if(.false.)),&
    !dir$ in(nnfirst: alloc_if(.true.) free_if(.false.)),&
    !dir$ in(nnlast: alloc_if(.true.) free_if(.false.)),&
    !dir$ in(q: alloc_if(.true.) free_if(.false.)),&
    !dir$ in(lj1: alloc_if(.true.) free_if(.false.)),&
    !dir$ in(lj2: alloc_if(.true.) free_if(.false.)),&
    !dir$ in(lj3: alloc_if(.true.) free_if(.false.)),&
    !dir$ in(lj4: alloc_if(.true.) free_if(.false.))


  end subroutine warm_up_mic_SP
  
  subroutine warm_up_mic_DP
    implicit none
    
    
    !dir$ offload_transfer target(mic:0) in(pdp: alloc_if(.true.) free_if(.false.)),&                                    
    !dir$ in(ffdp: alloc_if(.true.) free_if(.false.)),&                                                                                
    !dir$ in(nlist: alloc_if(.true.) free_if(.false.)),&                                                  
    !dir$ in(nnfirst: alloc_if(.true.) free_if(.false.)),&                                                                                         
    !dir$ in(nnlast: alloc_if(.true.) free_if(.false.)),&                                                                                              
    !dir$ in(qdp: alloc_if(.true.) free_if(.false.)),&                                                                                              
    !dir$ in(lj1dp: alloc_if(.true.) free_if(.false.)),&              
    !dir$ in(lj2dp: alloc_if(.true.) free_if(.false.)),&    
    !dir$ in(lj3dp: alloc_if(.true.) free_if(.false.)),&                                                                 
    !dir$ in(lj4dp: alloc_if(.true.) free_if(.false.))                                         

  end subroutine warm_up_mic_DP




  subroutine cool_down_mic_SP
    
    !dir$ offload_transfer target(mic:0) in(position: alloc_if(.false.) free_if(.true.)),&
    !dir$ in(ff: alloc_if(.false.) free_if(.true.)),&
    !dir$ in(nlist: alloc_if(.false.) free_if(.true.)),&
    !dir$ in(nnfirst: alloc_if(.false.) free_if(.true.)),&
    !dir$ in(nnlast: alloc_if(.false.) free_if(.true.)),&
    !dir$ in(q: alloc_if(.false.) free_if(.true.)),&
    !dir$ in(lj1: alloc_if(.false.) free_if(.true.)),&
    !dir$ in(lj2: alloc_if(.false.) free_if(.true.)),&
    !dir$ in(lj3: alloc_if(.false.) free_if(.true.)),&
    !dir$ in(lj4: alloc_if(.false.) free_if(.true.))
  end subroutine cool_down_mic_SP


  

  


 


 





  subroutine lj_cut_coul_dsf_newton_SP(step,printflag)
    implicit none
    real*4 :: force,forcelj,forcecoul
    real*4 :: x1,y1,z1,x2,y2,z2
    real*4 :: dx,dy,dz,dr,dr2,dr2i,dr6i,dr12i,r
    real*4 :: t,qtmp
    double precision :: ffx,ffy,ffz
    real*4 :: boxdx,boxdy,boxdz
    real*4 :: erfcc,erfcd
    real*4 :: prefactor
    real*4 :: coul_fudge, lj_fudge

    integer :: i,j,l,step
    integer :: itype,jtype,neigh
    integer :: tid,num
    integer :: offset
    integer :: T1,T2,clock_rate,clock_max
    integer :: printflag

    potential = 0.0d0; e_coul = 0.0d0


    
    
    call system_clock(T1,clock_rate,clock_max)

    !dir$ offload begin target(mic:0) in(position: alloc_if(.false.),free_if(.false.)),&
    !dir$ inout(ff: alloc_if(.false.), free_if(.false.)),&
    !dir$ in(nlist: alloc_if(.false.) free_if(.false.)),&
    !dir$ in(nnfirst: alloc_if(.false.) free_if(.false.)),&
    !dir$ in(nnlast: alloc_if(.false.) free_if(.false.)),&
    !dir$ in(q: alloc_if(.false.) free_if(.false.)),&
    !dir$ nocopy(lj1: alloc_if(.false.) free_if(.false.)),&
    !dir$ nocopy(lj2: alloc_if(.false.) free_if(.false.)),&
    !dir$ nocopy(lj3: alloc_if(.false.) free_if(.false.)),&
    !dir$ nocopy(lj4: alloc_if(.false.) free_if(.false.))

    !$omp parallel do schedule(dynamic) reduction(+:potential,e_coul,ffx,ffy,ffz) default(private),&
    !$omp& shared(box,ibox,rcut2,cut_coulsq,alpha,EWALD_P,MY_PIS,np,numAtomType),&
    !$omp& shared(A1,A2,A3,A4,A5),&
    !$omp& shared(e_shift,f_shift),&
    !$omp& shared(position,ff,nlist,nnfirst,nnlast,q,lj1,lj2,lj3,lj4)
    do i = 1 ,np
       x1 = position(i)%x; y1 = position(i)%y; z1 = position(i)%z; itype = position(i)%type
       qtmp = q(i)

       ffx = 0.0d0; ffy = 0.0d0; ffz = 0.0d0

   
       !dir$ simd reduction(+:potential,e_coul,ffx,ffy,ffz)
       do j= nnfirst(i),nnlast(i)
          
          neigh = nlist(j)
          x2 = position(neigh)%x; y2 = position(neigh)%y; z2 = position(neigh)%z; jtype = position(neigh)%type


          dx = x1-x2; dy = y1-y2; dz = z1-z2
         
          boxdx = dx*ibox; boxdy = dy*ibox; boxdz = dz*ibox
          boxdx = (boxdx+sign(1/(epsilon(boxdx)),boxdx)  )-sign(1/epsilon(boxdx),dx)
          boxdy = (boxdy+sign(1/(epsilon(boxdy)),boxdy)  )-sign(1/epsilon(boxdy),dy)
          boxdz = (boxdz+sign(1/(epsilon(boxdz)),boxdz)  )-sign(1/epsilon(boxdz),dz)

          dx = dx-box*boxdx; dy = dy-box*boxdy; dz = dz-box*boxdz

          dr2 = dx*dx + dy*dy + dz*dz

       
          coul_fudge = 0.0d0
          if(dr2.lt.cut_coulsq)coul_fudge = 1.0d0
                 
          r = sqrt(dr2)
          prefactor =  coul_fudge*qtmp*q(neigh)/r
          erfcd = exp(-alpha*alpha*r*r)
          t = 1.0 / (1.0 + EWALD_P*alpha*r)
          erfcc = t * (A1+t*(A2+t*(A3+t*(A4+t*A5)))) * erfcd
          forcecoul = prefactor * (erfcc/r + 2.0*alpha/MY_PIS * erfcd +&
               r*f_shift) * r

          e_coul = e_coul + prefactor*(erfcc-r*e_shift-dr2*f_shift)   

          
          lj_fudge = 0.0d0
          if(dr2.lt.rcut2)lj_fudge = 1.0d0
          
          dr2i = 1.0d0/dr2
          dr6i = dr2i*dr2i*dr2i
          
          offset  = numAtomType*(itype-1)+jtype
          forcelj = lj_fudge*dr6i*(lj1(offset)*dr6i-lj2(offset))
          
          potential = potential + lj_fudge*dr6i*(dr6i*lj3(offset)-lj4(offset))      
        

          force   = (forcecoul+forcelj)*dr2i
          
          ffx = ffx + dx*force
          ffy = ffy + dy*force
          ffz = ffz + dz*force
          
          
          ff(neigh)%x = ff(neigh)%x - dx*force
          ff(neigh)%y = ff(neigh)%y - dy*force
          ff(neigh)%z = ff(neigh)%z - dz*force   
       enddo

       ff(i)%x = ff(i)%x + ffx
       ff(i)%y = ff(i)%y + ffy
       ff(i)%z = ff(i)%z + ffz

    enddo
    !$omp end parallel do 
    !dir$ end offload

    call system_clock(T2,clock_rate,clock_max)


    if(printflag.eq.1)then
       print*,'elapsed time in DSF newton sp',real(T2-T1)/real(clock_rate)
       print*,'calculated potentials:',potential*reps, e_coul*reps
    endif

  end subroutine lj_cut_coul_dsf_newton_SP




  subroutine lj_cut_coul_dsf_newton_DP(step,printflag)
    implicit none
    double precision :: force,forcelj,forcecoul
    double precision :: x1,y1,z1,x2,y2,z2
    double precision :: dx,dy,dz,dr,dr2,dr2i,dr6i,dr12i,r
    double precision :: t,qtmp
    double precision :: ffx,ffy,ffz
    double precision :: boxdx,boxdy,boxdz
    double precision :: erfcc,erfcd
    double precision :: prefactor
    double precision :: coul_fudge,lj_fudge
    integer :: i,j,k,step
    integer :: itype,jtype, neigh
    integer :: tid,num
    integer :: offset
    integer :: T1,T2,clock_rate,clock_max
    integer :: printflag
    
    potential =0.0d0; e_coul = 0.0d0

    call system_clock(T1,clock_rate,clock_max)

  !  !dir$ offload begin target(mic:0) in(pdp: alloc_if(.false.),free_if(.false.)),&
  !  !dir$ inout(ffdp: alloc_if(.false.), free_if(.false.)),&
  !  !dir$ in(nlist: alloc_if(.false.) free_if(.false.)),&
  !  !dir$ in(nnfirst: alloc_if(.false.) free_if(.false.)),&
  !  !dir$ in(nnlast: alloc_if(.false.) free_if(.false.)),&
  !  !dir$ in(qdp: alloc_if(.false.) free_if(.false.)),&
  !  !dir$ nocopy(lj1dp: alloc_if(.false.) free_if(.false.)),&
  !  !dir$ nocopy(lj2dp: alloc_if(.false.) free_if(.false.)),&
  !  !dir$ nocopy(lj3dp: alloc_if(.false.) free_if(.false.)),&
  !  !dir$ nocopy(lj4dp: alloc_if(.false.) free_if(.false.))

  !  !$omp parallel do schedule(dynamic) reduction(+:potential,e_coul,ffx,ffy,ffz) default(private),&
  !  !$omp& shared(boxdp,iboxdp,rcut2,cut_coulsq,alphadp,EWALD_Pdp,MY_PISdp,np,numAtomType),&
  !  !$omp& shared(A1dp,A2dp,A3dp,A4dp,A5dp),&
  !  !$omp& shared(e_shiftdp,f_shiftdp),&
  !  !$omp& shared(pdp,ffdp,nlist,nnfirst,nnlast,qdp,lj1dp,lj2dp,lj3dp,lj4dp)
    do i = 1 ,np
       x1 = pdp(i)%x; y1 = pdp(i)%y; z1 = pdp(i)%z; itype = pdp(i)%type
       qtmp = qdp(i)

       ffx = 0.0d0; ffy = 0.0d0; ffz = 0.0d0

   
       !dir$ simd reduction(+:potential,e_coul,ffx,ffy,ffz)
       do j= nnfirst(i),nnlast(i)
          
          neigh = nlist(j)
          x2 = pdp(neigh)%x; y2 = pdp(neigh)%y; z2 = pdp(neigh)%z; jtype = pdp(neigh)%type


          dx = x1-x2; dy = y1-y2; dz = z1-z2
         
          boxdx = dx*iboxdp; boxdy = dy*iboxdp; boxdz = dz*iboxdp
          boxdx = (boxdx+sign(1/(epsilon(boxdx)),boxdx)  )-sign(1/epsilon(boxdx),dx)
          boxdy = (boxdy+sign(1/(epsilon(boxdy)),boxdy)  )-sign(1/epsilon(boxdy),dy)
          boxdz = (boxdz+sign(1/(epsilon(boxdz)),boxdz)  )-sign(1/epsilon(boxdz),dz)

          dx = dx-boxdp*boxdx; dy = dy-boxdp*boxdy; dz = dz-boxdp*boxdz

          dr2 = dx*dx + dy*dy + dz*dz

       
          coul_fudge = 0.0d0
          if(dr2.lt.cut_coulsq)coul_fudge = 1.0d0
                 
          r = sqrt(dr2)
          prefactor =  coul_fudge*qtmp*qdp(neigh)/r
          erfcd = exp(-alphadp*alphadp*r*r)
          t = 1.0 / (1.0 + EWALD_Pdp*alphadp*r)
          erfcc = t * (A1dp+t*(A2dp+t*(A3dp+t*(A4dp+t*A5dp)))) * erfcd
          forcecoul = prefactor * (erfcc/r + 2.0*alphadp/MY_PISdp * erfcd +&
               r*f_shiftdp) * r

          e_coul = e_coul + prefactor*(erfcc-r*e_shiftdp-dr2*f_shiftdp)   

          
          lj_fudge = 0.0d0
          if(dr2.lt.rcut2)lj_fudge = 1.0d0
          
          dr2i = 1.0d0/dr2
          dr6i = dr2i*dr2i*dr2i
          
          offset  = numAtomType*(itype-1)+jtype
          forcelj = lj_fudge*dr6i*(lj1dp(offset)*dr6i-lj2dp(offset))
          
          potential = potential + lj_fudge*dr6i*(dr6i*lj3dp(offset)-lj4dp(offset))      
          
          force   = (forcecoul+forcelj)*dr2i
          
          ffx = ffx + dx*force
          ffy = ffy + dy*force
          ffz = ffz + dz*force
          
          
          ffdp(neigh)%x = ffdp(neigh)%x - dx*force
          ffdp(neigh)%y = ffdp(neigh)%y - dy*force
          ffdp(neigh)%z = ffdp(neigh)%z - dz*force   
       enddo

       ffdp(i)%x = ffdp(i)%x + ffx
       ffdp(i)%y = ffdp(i)%y + ffy
       ffdp(i)%z = ffdp(i)%z + ffz

    enddo
 !   !$omp end parallel do 
 !   !dir$ end offload

    call system_clock(T2,clock_rate,clock_max)

    if(printflag.eq.1)then
       print*,'elapsed time in DSF newton dp',real(T2-T1)/real(clock_rate)
       print*,'calculated potentials:',potential*reps, e_coul*reps
    endif

  end subroutine lj_cut_coul_dsf_newton_DP


end module mod_force
