       Program cluster 
!****Solves classical equations of motion by symplectic integration*****
!*******Splits hamiltonian on Coulomb pairs and pairs in the laser******
!**************and the rest independent on kinetic energy***************
       integer outel(512)
       integer eiasm(512)
       real*8 Pi,co,dist
       real*8 tlen
       real*8  acci(64,3)
       real*8 acce(512,3)
       real*8 veli(64,3)
       real*8 vele(512,3)
       real*8 xyzi(64,3)
       real*8 xyze(512,3)
       real*8  eeintfu(512,512)
       real*8 iiintfu(64,64)
       real*8 iieintfu(64,512)
       real*8 eiintfu(512,64)
       real*8  dt
       integer  ionhist(512)
       real*8 ionchar(64)
	   integer epairlist(512,2)
	   integer ifpair
	   integer nepair
	   integer unpair
	   real*8 tempt
	   real*8 ifeiint0(512,64)
	   real*8 omega
	   real*8 E0(3)
           real*8 E00(3)
	   real*8 sins
	   real*8 tt
	   real*8 ww
	   real*8 corechar
	   real*8 mu
	   real*8 mi 
	   real*8 me 
	   real*8 dt0
           real*8 sumionchar
           real*8 x(900000),y(900000),z(900000)
           real*8 xi(80000),yi(80000),zi(80000) 
           integer nele
           integer nion
           integer nelat
           integer ntime
           integer ifcol(512,64)
           real*8 eee(16) 
           real*8 kkk(16)
           real*8 ener(16)

           real*8 xr(10,64)
           real*8 xrt(10)
           real*8 xrts(10)

cc       write(1,*) 200 ,0, 200
cc       write(1,*) 3
cc       write(1,*) 1.0
cc       do i=1,200
cc       write(1,*) 0.1 
cc       enddo
       Pi=3.14159265358979311
       tempt=0.0   
       nele=1
       nion=2
!       do i=1,512
       do i=1,nele
       ionhist(i)=1
       enddo
!          do i=1,64
           do i=1,nion
	   ionchar(i)=1.0
	   enddo
	   corechar=8.0
	   mi=20.1830*1836.0
	   me=1.0
	   omega=1.0/17.123
	   ww=omega
           dt0=3.141*2.0/omega
           tlen=2*3.141/omega*10.0
           dt=dt0
!*******************sets the initial conditions ***********************
!*****defines electron assigment to ions and time average initial******
!************interations in a classical hartree-fock way***************
!            do i=1,512
           do i=1,nele
!           do j=1,64
           do j=1,nion
	   ifeiint0(i,j)=1.0
           ifcol(i,j)=0
	   enddo
           eiasm(i)=1
	   enddo

!          do i=1,512
           do i=1,nele
       vele(i,1)=0.0
       vele(i,2)=0.0
       vele(i,3)=0.0
       xyze(i,1)=0.0
       xyze(i,2)=0.0
       xyze(i,3)=0.0
       enddo
       mu=mi
!      do i=1,64
       do i=1,nion
       veli(i,1)=0.0
       veli(i,2)=0.0
       veli(i,3)=0.0
       xyzi(i,1)=0.0
       xyzi(i,2)=0.0
       xyzi(i,3)=0.0
       enddo
cc      do i=1,3
cc      E00(i)=1.0
cc      enddo
       E00(1)=0.37*100.0
cc       E00(1)=0.0
       E00(2)=0.0
       E00(3)=0.0
       do i=1,1 
       vele(i,1)=-7.0
       vele(i,2)=-0.2
       vele(i,3)=-0.9
       xyze(i,1)=0.10
       xyze(i,2)=0.100000000000001
       xyze(i,3)=0.0
       enddo
       
       do i=2,2
       vele(i,1)=-0.8
       vele(i,2)=-1.2
       vele(i,3)=-0.2
       xyze(i,1)=4.00
       xyze(i,2)=2.000000000000001
       xyze(i,3)=0.0
       enddo


 
       do i=1,nion
       veli(i,1)=0.0
       veli(i,2)=0.0
       veli(i,3)=0.0
       xyzi(i,1)=0.0
       xyzi(i,2)=0.0
       xyzi(i,3)=0.0
       enddo
       eiasm(1)=1

cc       do i=1,3
cc       do j=1,3
cc       do k=1,3
cc       l=(i-1)*9+(j-1)*3+k
cc       xyzi(l,1)=-1.5*6.0+Qfloat(i-1)*6.0
cc       xyzi(l,2)=-1.5*6.0+Qfloat(j-1)*6.0
cc       xyzi(l,3)=-1.5*6.0+Qfloat(k-1)*6.0
cc       write(*,*) l
cc       write(*,*) xyzi(l,1)
cc       write(*,*) xyzi(l,2)
cc       write(*,*) xyzi(l,3)
cc       enddo
cc       enddo
cc       enddo
cc       stop
ccc       xyzi(1,1)=0.0
ccc       xyzi(1,2)=0.0
ccc       xyzi(1,3)=0.0
ccc       veli(1,1)=0.0

cc       xyzi(2,1)=1.0
cc       xyzi(2,2)=0.0
cc       xyzi(2,3)=0.0

       ntime=1
       nelat=8
       nion=27
       open(2,file='fort.2',status='old')
       do i=1,nion
       read(2,20) xyzi(i,1),xyzi(i,2),xyzi(i,3)
       enddo
20     format(3F15.8)
       close(2)

cc       xyzi(1,1)=0.0
cc       xyzi(1,2)=0.0
cc       xyzi(1,3)=0.0
 
       do i=1,nelat
       eee(i)=0.05
       kkk(i)=float(i)
       ener(i)=-float(i)**2.0/(2.0*1.0**2.0)
       write(*,*) ener(i)
       enddo  
       call prepfull(eee,ener,kkk,nelat,nion,xyze,vele,veli,xyzi,
     & ionchar,ifeiint0,eiasm,outel,ionhist)
       nele=nion*nelat
cc       nele=2
cc       nion=1
cc       do i=1,nele
cc       write(2,20) xyze(i,1),xyze(i,2),xyze(i,3)
cc       write(2,*) eiasm(i),outel(i),ionhist(i)
cc20     FORMAT(3F15.8)
cc       enddo 
cc       open(5,file='plot.dat')
cc       do i=1,nion
cc       write(5,21) xyzi(i,1),xyzi(i,2),xyzi(i,3)
cc21     FORMAT(3F15.8)
cc       enddo


cc       do i=1,nele
cc       write(5,30) xyze(i,1),xyze(i,2),xyze(i,3)
cc30     FORMAT(3F15.8)
cc       enddo
cc       close(5)
cc       stop
!************************main time-loop here****************************
       tt=0.0
       
       open(7,file='plot2.dat')
       open(8,file='plot3.dat')
       open(9,file='plot4.dat')
      n=0
      xrts(1)=0.0
      xrts(2)=0.0
10    continue
cc       do n=1,ntime
       do i=1,3
cc       E0(i)=E00(i)*EXP(-(tt-tlen/2.0)**2.0/(tlen/2.4)**2.0)
       E0(i)=E00(i)*(SIN(tt/tlen*2.0*Pi-0.5*Pi)+1.0)/2.0
       enddo
!10     continue
       write(*,*) n,'hi'
       ifpair=1
       write(*,*) E0(1),'E0'
!**********************N-N Coulomb Drift C(t/2)*************************
!*********************pair assigment only here !!!**********************
       Call coulombmain(xyze,xyzi,veli,vele,ionhist,ionchar,
     & epairlist,dt,ifpair,nepair,unpair,tempt, ifeiint0,omega,E0,eiasm,
     & nele, nion,tt,outel)
cc       write(*,*) dt, 'after 1st coulomb'
cc       write(*,*) nepair,'nepair',unpair
cc       do i=1,nepair
cc       write(*,*) epairlist(i,1),epairlist(i,2)
cc       enddo 
cc       do i=1,nion
cc       write(*,*) ionchar(i)
cc       enddo
cc       do i=i,nele
cc       write(*,*) eiasm(i)
cc       enddo
cc       stop
cc       write(*,*) dt,'dt'
cc       write(*,*) ifeiint0(1,1), 'why am i changing',eiasm(1)
!********************residual Kick K(t)*********************************
!*********all particles interact with the laser*************************
!*ions are coupled by charge of a nucleus since electrons individually**
!**otherwise should be on average only which would prevent ionization***
!***********************************************************************
       
cc       do i=1,nele
cc       write(*,*) xyze(i,1),xyze(i,2),xyze(i,3)
cc       enddo
cc       do i=1,nion
cc       write(*,*) xyzi(i,1),xyzi(i,2),xyzi(i,3)
cc       enddo
cc       do i=1,nele
cc       write(*,*) xyze(i,1),xyze(i,2),xyze(i,3)
cc       enddo
       
cc       do i=1,nele
cc       write(*,*) eiasm(i)
cc       enddo
cc       stop

cc        write(3,*) 'before  acce'
cc        do i=1,nion
cc        write(3,*) xyzi(i,1),xyzi(i,2),xyzi(i,3)
cc        enddo
cc        do i=1,nele
cc        write(3,*) xyze(i,1),xyze(i,2),xyze(i,3)
cc        enddo
 
         call  accels(xyze,xyzi,acce,acci,ionhist,ionchar,epairlist,
     &          ifeiint0,nepair,eiasm,nele,nion,outel,corechar)
cc        write(3,*) 'after acce'
cc        do i=1,nion
cc        write(3,*) acci(i,1),acci(i,2),acci(i,3)
cc        enddo
cc        do i=1,nele
cc        write(3,*) acce(i,1),acce(i,2),acce(i,3)
cc        enddo
cc         write(3,*) 'after  full'
cc        do i=1,nion
cc        write(3,*) veli(i,1),veli(i,2),veli(i,3)
cc        enddo
cc        do i=1,nele
cc        write(3,*) vele(i,1),vele(i,2),vele(i,3)
cc        enddo
 


cc       write(*,*) 'aciiiiiiiiiiiiiiiiiiiiiiiiii'

       sins=Sin(ww*tt)-Sin(ww*(tt+2.0*dt))
cc       sins=ww*2.0*dt
       mu=me
!       do i=1,512
       do i=1,nele
       co=1.0
cc       if(ionhist(i).eq.0) then
cc       dist=ABS((xyzi(eiasm(i),1)-xyze(i,1))
cc     & **2.0+
cc     & (xyzi(eiasm(i),2)-xyze(i,2))**2.0+(xyzi(eiasm(i),3)-
cc     & xyze(i,3))**2.0)**(1.0/2.0)
cc       if(dist.lt.0.1) co=0.0 
cc       endif
       vele(i,1)=vele(i,1)+2.0*dt*acce(i,1)+(E0(1)/mu/ww)*sins*co
       vele(i,2)=vele(i,2)+2.0*dt*acce(i,2)+(E0(2)/mu/ww)*sins*co
       vele(i,3)=vele(i,3)+2.0*dt*acce(i,3)+(E0(3)/mu/ww)*sins*co
       enddo
       mu=mi
!      do i=1,64
       do i=1,nion
       veli(i,1)=veli(i,1)+2.0*dt*acci(i,1)-(E0(1)/mu/ww)*sins*corechar
       veli(i,2)=veli(i,2)+2.0*dt*acci(i,2)-(E0(2)/mu/ww)*sins*corechar
       veli(i,3)=veli(i,3)+2.0*dt*acci(i,3)-(E0(3)/mu/ww)*sins*corechar
       enddo
cc       write(3,*) 'after  full'
cc        do i=1,nion
cc        write(3,*) veli(i,1),veli(i,2),veli(i,3)
cc        enddo
cc        do i=1,nele
cc        write(3,*) vele(i,1),vele(i,2),vele(i,3)
cc        enddo

!****correction since electron pairs were in the laser field already****
!****************and  e-e scattering was  solved exactly****************
       mu=me
       do i=1,nepair
       if((epairlist(i,2).lt.1000).and.(epairlist(i,1).lt.1000)) then
       vele(epairlist(i,1),1)=vele(epairlist(i,1),1)-(E0(1)/mu/ww)*sins
       vele(epairlist(i,1),2)=vele(epairlist(i,1),2)-(E0(2)/mu/ww)*sins
       vele(epairlist(i,1),3)=vele(epairlist(i,1),3)-(E0(3)/mu/ww)*sins
       vele(epairlist(i,2),1)=vele(epairlist(i,2),1)-(E0(1)/mu/ww)*sins
       vele(epairlist(i,2),2)=vele(epairlist(i,2),2)-(E0(2)/mu/ww)*sins
       vele(epairlist(i,2),3)=vele(epairlist(i,2),3)-(E0(3)/mu/ww)*sins
       endif
       enddo
!********************or to unpaired particle if any*********************
cc       unpair=0
       if (unpair.gt.0) then
	   i=unpair
       if (unpair.gt.1000) then   
       mu=mi
       veli(i-1000,1)=veli(i-1000,1)+(E0(1)/mu/ww)*sins*corechar
       veli(i-1000,2)=veli(i-1000,2)+(E0(2)/mu/ww)*sins*corechar
       veli(i-1000,3)=veli(i-1000,3)+(E0(3)/mu/ww)*sins*corechar
       
       else  
       mu=me
       vele(i,1)=vele(i,1)-(E0(1)/mu/ww)*sins
       vele(i,2)=vele(i,2)-(E0(2)/mu/ww)*sins
       vele(i,3)=vele(i,3)-(E0(3)/mu/ww)*sins
       
       endif    
         
       endif
      
 

!**********************N-N Coulomb Drift C(t/2)*************************
       ifpair=0
cc       write(*,*) 'before calling',unpair
cc       write(*,*) dt , 'before 2nd coulomb'
       Call coulombmain(xyze,xyzi,veli,vele,ionhist,ionchar,
     & epairlist,dt,ifpair,nepair,unpair,tempt,ifeiint0,omega,E0,eiasm,
     & nele,nion,tt,outel)
 

cc        write(3,*) 'after  full'
cc        do i=1,nion
cc        write(3,*) xyzi(i,1),xyzi(i,2),xyzi(i,3)
cc        enddo
cc        do i=1,nele
cc        write(3,*) xyze(i,1),xyze(i,2),xyze(i,3)
cc        enddo

cc       write(*,*) dt , 'after 2nd coulomb'
       write(*,*) 'after 2nd' 
       tt=tt+dt*2.0
cc       write(*,*) tt/dt0
       write(*,*) tt*omega/(2.0*3.141)
cc     if (tt.gt.tlen) stop
cc       goto 10 
c       write(1,*) xyze(1,1)
c       write(1,*) xyze(1,2) 
c       write(1,*) xyze(1,3)
cc       do k=1,nele
cc       m=(n-1)*nele+k
cc       x(m)= xyze(k,1)
cc       y(m)= xyze(k,2)
cc       z(m)= xyze(k,3)
cc       enddo
cc       do k=1,nion
cc       m=(n-1)*nion+k
cc       xi(m)= xyzi(k,1)
cc       yi(m)= xyzi(k,2)
cc       zi(m)= xyzi(k,3)
cc       enddo
       call   xray(xyze,vele,veli,xyzi,xr,xrt,ifcol,nele,nion,
     &    ionhist)
       xrts(1)=xrts(1)+xrt(1)
       xrts(2)=xrts(2)+xrt(2)
       write(*,*) sins,'*********%%%%%%%%%%%%%*******'
       IF(ABS((n-1)/300-(n-1)/300.).LT.1.E-11) THEN
       do i=1,nion
       write(7,23) xyzi(i,1),xyzi(i,2),xyzi(i,3),ionchar(i)*1.0
23     FORMAT(4F15.8)
       enddo
       sumionchar=0.0 
       do i=1,nion
       sumionchar=sumionchar+ionchar(i)-1.0
       enddo
       write(8,24) tt,sumionchar/float(nion),E0(1)
24     FORMAT(3F15.8)
cc       write(7,24) tt,ionchar(1)
cc24     FORMAT(2F15.8)
cc       do i=1,nele
       i=4
       write(7,34) xyze(i,1),xyze(i,2),xyze(i,3),
     & (vele(i,1)**2.0+vele(i,2)**2.0+vele(i,3)**2.0)/2.0
       write(9,35) tt,E0(1)/mu/ww*sins,
     & (xyze(1,1)**2.0+xyze(1,2)**2.0+xyze(1,3)**2.0)**0.5,
     & (xyze(2,1)**2.0+xyze(2,2)**2.0+xyze(2,3)**2.0)**0.5,
     & (xyze(3,1)**2.0+xyze(3,2)**2.0+xyze(3,3)**2.0)**0.5,
     & (xyze(4,1)**2.0+xyze(4,2)**2.0+xyze(4,3)**2.0)**0.5,
     & (xyze(5,1)**2.0+xyze(5,2)**2.0+xyze(5,3)**2.0)**0.5,
     & (xyze(6,1)**2.0+xyze(6,2)**2.0+xyze(6,3)**2.0)**0.5,
     & (xyze(7,1)**2.0+xyze(7,2)**2.0+xyze(7,3)**2.0)**0.5,
     & (xyze(8,1)**2.0+xyze(8,2)**2.0+xyze(8,3)**2.0)**0.5
cc     & (vele(1,1)**2.0+vele(1,2)**2.0+vele(1,3)**2.0)**0.5,
cc     & (vele(2,1)**2.0+vele(2,2)**2.0+vele(2,3)**2.0)**0.5,
cc     & (vele(3,1)**2.0+vele(3,2)**2.0+vele(3,3)**2.0)**0.5,
cc     & (vele(4,1)**2.0+vele(4,2)**2.0+vele(4,3)**2.0)**0.5,
cc     & (vele(5,1)**2.0+vele(5,2)**2.0+vele(5,3)**2.0)**0.5,
cc     & (vele(6,1)**2.0+vele(6,2)**2.0+vele(6,3)**2.0)**0.5,
cc     & (vele(7,1)**2.0+vele(7,2)**2.0+vele(7,3)**2.0)**0.5,
cc     & (vele(8,1)**2.0+vele(8,2)**2.0+vele(8,3)**2.0)**0.5

cc     &  (vele(1,1)**2.0+vele(1,2)**2.0+vele(1,3)**2.0)/2.0+
cc     & -1.0*(xyze(1,1)**2.0+xyze(1,2)**2.0+xyze(1,3)**2.0)**(-0.5),
cc     &  (vele(2,1)**2.0+vele(2,2)**2.0+vele(2,3)**2.0)/2.0+
cc     & -2.0*(xyze(2,1)**2.0+xyze(2,2)**2.0+xyze(2,3)**2.0)**(-0.5),
cc     &  (vele(3,1)**2.0+vele(3,2)**2.0+vele(3,3)**2.0)/2.0+
cc     & -3.0*(xyze(3,1)**2.0+xyze(3,2)**2.0+xyze(3,3)**2.0)**(-0.5),
cc     &  (vele(4,1)**2.0+vele(4,2)**2.0+vele(4,3)**2.0)/2.0+
cc     & -4.0*(xyze(4,1)**2.0+xyze(4,2)**2.0+xyze(4,3)**2.0)**(-0.5),
cc     &  (vele(5,1)**2.0+vele(5,2)**2.0+vele(5,3)**2.0)/2.0+
cc     & -5.0*(xyze(5,1)**2.0+xyze(5,2)**2.0+xyze(5,3)**2.0)**(-0.5),
cc     &  (vele(6,1)**2.0+vele(6,2)**2.0+vele(6,3)**2.0)/2.0+
cc     & -6.0*(xyze(6,1)**2.0+xyze(6,2)**2.0+xyze(6,3)**2.0)**(-0.5),
cc     &  (vele(7,1)**2.0+vele(7,2)**2.0+vele(7,3)**2.0)/2.0+
cc     & -7.0*(xyze(7,1)**2.0+xyze(7,2)**2.0+xyze(7,3)**2.0)**(-0.5),
cc     &  (vele(8,1)**2.0+vele(8,2)**2.0+vele(8,3)**2.0)/2.0+
cc     & -8.0*(xyze(8,1)**2.0+xyze(8,2)**2.0+xyze(8,3)**2.0)**(-0.5)



34     FORMAT(4F15.8)
35     FORMAT(10F15.8)
cc       enddo

       ENDIF 

cc       enddo
       if (tt.gt.tlen) stop
       n=n+1 
       write(*,*) n,'nn'
       write(*,29) xrts(1),xrts(2)
29     FORMAT(2F30.15)
       goto 10
!************main loop ends here****************************************
       nele=1
       nion=1
       write(1,*) (nele+nion)*ntime ,nion*ntime, nele*ntime
       write(1,*) 3
       write(1,*) 1.0
       do i=1,(nele+nion)*ntime
       write(1,*) 0.1
       enddo

       do i=1,nion*ntime
       write(1,*) xi(i)
       enddo
       do i=1,nele*ntime
       write(1,*) x(i)
       enddo
       
       do i=1,nion*ntime
       write(1,*) yi(i)
       enddo
       do i=1,nele*ntime
       write(1,*) y(i)
       enddo  

       do i=1,nion*ntime
       write(1,*) zi(i)
       enddo
       do i=1,nele*ntime
       write(1,*) z(i)
       enddo  
     
       do i=1,nion*ntime
       write(1,*) xi(i)
       enddo
       do i=1,nele*ntime
       write(1,*) x(i)
       enddo

       do i=1,nion*ntime
       write(1,*) yi(i)
       enddo
       do i=1,nele*ntime
       write(1,*) y(i)
       enddo
 
       do i=1,nion*ntime
       write(1,*) zi(i)
       enddo
       do i=1,nele*ntime
       write(1,*) z(i)
       enddo




       do i=1,nele*ntime
       write(1,*) 0.1
       enddo
       
       do i=1,4*nion*ntime
       write(1,*) 0.1
       enddo

       do i=1,nele*ntime
       write(1,*) 0.1
       enddo
       do i=1,nele*ntime
       write(1,*) 0.0
       enddo
       do i=1,(nele+nion)*ntime
       write(1,*) 0.1 
       enddo




       

!***********************************************************************  
      
       open(6,file='plot1.dat')
       do i=1,nion
       write(6,22) xyzi(i,1),xyzi(i,2),xyzi(i,3)
22     FORMAT(3F15.8)
       enddo


       do i=1,nele
       write(6,33) xyze(i,1),xyze(i,2),xyze(i,3)
33     FORMAT(3F15.8)
       enddo
       close(6)
       close(7)
       close(8)
       close(9)
       end 

       subroutine  coulombmain(xyze,xyzi,veli,vele,ionhist,ionchar
     & ,epairlist,dt,ifpair,nepair,unpair,tempt, 
     &  ifeiint0,omega,E0,eiasm,nele,nion,tt,outel)
!******************************************************************************** 
!**solves the Coulomb problem exactly for strongest interacting pairs************
!***************pairing should occure only once per full split operator !!!******
!*************and therefore determination of ionization conditions***************
       integer outel(512)
       integer nele
       integer nion
       real*8 acci(64,3)
       real*8 acce(512,3)
       real*8 veli(64,3)
       real*8 velio(64,3)
       real*8 vele(512,3)
       real*8 xyzi(64,3)
       real*8 xyzio(64,3) 
       real*8 xyze(512,3)
       real*8 eeint(512,512)
       real*8 iiint(64,64)
       real*8 ieint(64,512)
       real*8 eiint(512,64)
       real*8 eeintfu(512,512)
       real*8 iiintfu(64,64)
       real*8 ieintfu(64,512)
       real*8 eiintfu(512,64)
       real*8 eedist(512,512)
       real*8 iidist(64,64)
       real*8 iedist(64,512)
       real*8 eidist(512,64)
       real*8 ifeeint(512,512)
       real*8 ifiiint(64,64)
       real*8 ifieint(64,512)
       real*8 ifeiint(512,64)
       real*8 ifeiint0(512,64)
       real*8 atomsize
       real*8 tempee(512,512)
       real*8 tempei(512,64)
       real*8 tempie(64,512)
       real*8 tempii(64,64)
       real*8 E0(3)
       real*8 omega
       integer ionhist(512)
       integer ionhist1(512)
       integer eiasm(512)
       integer epairlist(512,2)
       integer ipairlist(512,2)
       integer              h
       real*8  ionchar(64)
       real*8      ioncharm(64,512)
       real*8                        zion
       real*8                        dt
       real*8                        mu
       integer                     ninte
	   integer                     ninti
	   integer                     ifpair	   
       integer     intel(512)
	   integer     intil(64)
	   integer     nepair
	   integer     zn1
	   real*8        aa 
	   integer     zn
	   real*8        dist 
	   integer     unpair 
	   real*8        mindt
	   real*8        dt0
	   real*8        dist2
	   real*8        tempdt
	   real*8        nstep
	   real*8        nmax
	   real*8        tempt
	   real*8        xx
	   real*8        yy
	   real*8        zz
	   real*8        px
	   real*8        py
	   real*8        pz
	   real*8        me
	   real*8        mi
	   real*8        m1
	   real*8        m2
	   real*8        kk
	   real*8        Rxx
	   real*8        Ryy
	   real*8        Rzz
	   real*8        Vxx
	   real*8        Vyy
	   real*8        Vzz
	   real*8        ww
           real*8        tt	   
           real*8        E00(3)
           real*8        corechar 
           corechar=8.0
           atomsize=5.67
           E00(1)=0.0
           E00(2)=0.0
	   E00(3)=0.0 
           me=1.0
	   mi=20.1830*1836.0
	   ww=omega
           dt0=2*3.141/omega/64.0/17.0
           nstep=10.0
           ninte=0
       if(ifpair.eq.1) then 



!**************Ionization determination*********************************
!*****************Moved from accels*************************************
!**************if ionized***********************************************
!******************distance from assigned ions**************************
       do  i=1,nion
       do  j=1,nele
       if(i.eq.eiasm(j)) then 
cc       if((i.eq.eiasm(j)).and.(ionhist(j).eq.0.0)) then
       eidist(j,i)=ABS((xyzi(i,1)-xyze(j,1))**
     &2.0+(xyzi(i,2)-xyze(j,2))**2.0+(xyzi(i,3)-xyze(j,3))
     &**2.0)**(1.0/2.0)
       endif
       enddo
       enddo

       do  j=1,nele
       if((ionhist(j).eq.0).and.(eidist(j,eiasm(j)).gt.atomsize))
     & ionhist(j)=1
       
       if((ionhist(j).eq.2)) ionhist(j)=1

       if(eidist(j,eiasm(j)).gt.3000.0) ionhist(j)=2
       enddo
!**************Charge of each ion***************************************
!       do  i=1,64
       do  i=1,nion
       ionchar(i)=1.0
!      do  j=1,512
       do  j=1,nele
       if((i.eq.eiasm(j)).and.(ionhist(j).gt.0))
     &ionchar(i)=1.0+ionchar(i)
cc       write(*,*) nele,ionhist(j)
cc       write(*,*) j, eiasm(j)
       enddo
       if(ionchar(i).gt.corechar) ionchar(i)=corechar
cc       write(*,*) 'in accels'
cc       write(*,*) ionchar(i)
       enddo


!**********************Outer electrons**********************************
       do i=1,nele
       outel(i)=0
cc       write(*,*) ifeiint0(i,eiasm(i)),ionchar(eiasm(i)),
cc     &  eiasm(i),i,'**********'
       if((ifeiint0(i,eiasm(i)).eq.ionchar(eiasm(i))).and.
     &   (ionhist(i).eq.0))  outel(i)=1
       enddo


cc       write(*,*) ifpair , 'ifpair' 
       j=0
!       do i=1,512
       do i=1,nele
cc       write(*,*) ionhist(i),i,outel(i)
cc       write(*,*) eiasm(i)
       if (ionhist(i).eq.1) then
       j=j+1
       intel(j)=i
       endif
       if ((ionhist(i).eq.0).and.(outel(i).eq.1)) then
       j=j+1
       intel(j)=i
       endif
cc       eiasm(i)=1
       enddo
       ninte=j
cc       write(*,*) ninte,nion ,nele,ninti,'aaa'
       write(*,*) ninte,nepair
       j=0
!      do i=1,64
       do i=1,nion
       if (ionchar(i).gt.0.0) then
cc       write(*,*) ionchar(i),'ionchaaaaaaaaaaaaaaaaaaaaaaaaaaaaaar'
       j=j+1
       intil(j)=i
       endif
       enddo
       ninti=j
cc       write(*,*) 'after calculating',ninti
cc       write(3,*) 'next before pairing'  
cc       do i=1,ninti
cc       write(3,*) xyzi(intil(i),1),xyzi(intil(i),2),xyzi(intil(i),3) 
cc       enddo 
cc       write(3,*) 'elec'
cc       do i=1,ninte
cc       write(3,*) xyze(intel(i),1),xyze(intel(i),2),xyze(intel(i),3)
cc       enddo
!******************************Pairing particles*********************************
!********according to closest "charge distance" and natural order****************
!*****************************this is kind of random*****************************
!*******probably should be replaced by quick sorting of all dists****************
!*********************and second linear sorting of distinct pairs****************
!**************extra inner loop versus all distance calculation******************
!***********************Electrons and ions **************************************
       nepair=0
!********************************************************************************
!**************************electron first in pair********************************
!********************************************************************************
       do i=1,ninte
       zn1=0
       aa=10000000000.0
!*************************looking around electrons*******************************
       do j=1,ninte
       zn=0
       do k=1,nepair
       if((intel(j).eq.epairlist(k,1)).or.
     & (intel(j).eq.epairlist(k,2))) zn=1
       if((intel(i).eq.epairlist(k,1)).or.
     & (intel(i).eq.epairlist(k,2))) zn=1
       enddo
!**********************outer cannot be paired************************************
       if(ionhist(intel(j)).eq.0) zn=1
       if(ionhist(intel(i)).eq.0) zn=1
c       write(*,*) 'checking zn', zn 
       if (zn.eq.0) then 
       dist=ABS((xyze(intel(i),1)-xyze(intel(j),1))
     & **2.0+
     & (xyze(intel(i),2)-xyze(intel(j),2))**2.0+(xyze(intel(i),3)-
     & xyze(intel(j),3))**2.0)**(1.0/2.0)
       if((aa.gt.dist).and.(i.ne.j)) then
       epairlist(nepair+1,1)=intel(i)
       epairlist(nepair+1,2)=intel(j)
       aa=dist
       zn1=1
       endif
       endif
c       write(*,*) 'in e-e loop',zn1
       enddo
c       write(*,*) 'nepair after e-e loop', nepair
!************************looking around ions ***********************************
       do j=1,ninti
       zn=0
       do k=1,nepair
       if((intil(j).eq.epairlist(k,1)-1000).or.
     & (intil(j).eq.epairlist(k,2)-1000))  zn=1
       if((intel(i).eq.epairlist(k,1)).or.
     & (intel(i).eq.epairlist(k,2)))  zn=1
       enddo
!**********************outer cannot be paired************************************
       if(ionhist(intel(i)).eq.0) zn=1

c      write(*,*) 'why did i change after loop', nepair
c      write(*,*) 'checking zn', zn 
       if (zn.eq.0) then 
       dist=ABS((xyze(intel(i),1)-xyzi(intil(j),1))
     & **2.0+
     & (xyze(intel(i),2)-xyzi(intil(j),2))**2.0+(xyze(intel(i),3)-
     & xyzi(intil(j),3))**2.0)**(1.0/2.0)/(ionchar(j)**0.5)
       if(aa.gt.dist) then
       epairlist(nepair+1,1)=intel(i)
       epairlist(nepair+1,2)=intil(j)+1000
       aa=dist
       zn1=1
c       write(*,*) 'why did i change', nepair
       endif
       endif
      
c       write(*,*) 'in e-i loop',zn1,nepair
       enddo
c       write(*,*) 'number of pairs before',nepair
        
       if(zn1.eq.1) nepair=nepair+1
       enddo
cc       write(7,*) 'number of paires first', nepair
!********************************************************************************
!*******************************ion first in pair********************************
!********************************************************************************
       do i=1,ninti
cc       write(*,*) i
       zn1=0
       aa=10000000000.0
!*************************looking around electrons*******************************
       do j=1,ninte
cc       write(*,*) 'am i in electrons'
       zn=0
       do k=1,nepair
       if((intel(j).eq.epairlist(k,1)).or.
     & (intel(j).eq.epairlist(k,2))) zn=1
       if((intil(i).eq.epairlist(k,1)-1000).or.
     & (intil(i).eq.epairlist(k,2)-1000)) zn=1
       enddo
!**********************outer cannot be paired************************************
       if(ionhist(intel(j)).eq.0) zn=1
       
       if (zn.eq.0) then 
       dist=ABS((xyzi(intil(i),1)-xyze(intel(j),1))
     & **2.0+
     & (xyzi(intil(i),2)-xyze(intel(j),2))**2.0+(xyzi(intil(i),3)-
     & xyze(intel(j),3))**2.0)**(1.0/2.0)/(ionchar(i)**0.5)
       if(aa.gt.dist) then
       epairlist(nepair+1,1)=intil(i)+1000
       epairlist(nepair+1,2)=intel(j)
       aa=dist
       zn1=1
       endif
       endif
cc       write(3,*) 'in i-e loop',zn1,dist
       enddo
!************************looking around ions ************************************
       do j=1,ninti
       if(ionchar(j).ne.0) then
       zn=0
       do k=1,nepair
       if((intil(j).eq.epairlist(k,1)-1000).or.
     & (intil(j).eq.epairlist(k,2)-1000)) zn=1
       if((intil(i).eq.epairlist(k,1)-1000).or.
     & (intil(i).eq.epairlist(k,2)-1000)) zn=1
       enddo
       
       if (zn.eq.0) then 
       dist=ABS((xyzi(intil(i),1)-xyzi(intil(j),1))
     & **2.0+
     & (xyzi(intil(i),2)-xyzi(intil(j),2))**2.0+(xyzi(intil(i),3)-
     & xyzi(intil(j),3))**2.0)**(1.0/2.0)/
     &(ionchar(i)**0.5*ionchar(j)**0.5)
cc       write(*,*) ionchar(i)**0.5*ionchar(j)**0.5, 'ion'
       if((aa.gt.dist).and.(i.ne.j)) then
       epairlist(nepair+1,1)=intil(i)+1000
       epairlist(nepair+1,2)=intil(j)+1000
cc       write(*,*) epairlist(nepair+1,1),epairlist(nepair+1,2),'opop'
       aa=dist
       zn1=1
       endif
       endif
       endif
c       write(*,*) 'in i-i loop',zn1 
       enddo
       
       
       if(zn1.eq.1) nepair=nepair+1
ccc       write(*,*) nepair
       enddo
cc       do i=1,nepair
cc       write(4,*) epairlist(i,1),epairlist(i,2),i
cc       enddo
cc       write(*,*) 'number of pairs', nepair
!********************looks for unpaired interacting particle*********************
!*************************checks pairing*****************************************
!********all interacting particles must be paired or except one******************
       if((nepair.ne.(ninte+ninti)/2).and.(nepair.ne.(ninte+ninti-1)/2)) 
     & write(*,*) 'pairing failed',ninte,ninti,nepair
       unpair=0  
!***********extracts left unpaired interacting particle if any*******************
cc       if(nepair.eq.(ninte+ninti-1)/2) then 
       write(*,*) 'am i checking'
!****************************if unpaired is electron*****************************
       do i=1,ninte
       zn=0
       do k=1,nepair
       if((intel(i).eq.epairlist(k,1)).or.(intel(i).eq.epairlist(k,2))) 
     & zn=1
cc       write(*,*) intel(i), epairlist(k,1), epairlist(k,2)
       enddo
       if (ionhist(intel(i)).eq.0) zn=1
       if(zn.eq.0) unpair=intel(i)
       enddo
       write(*,*) 'after unpared is electron' , unpair,ninte,nepair
!*******************************if unpaired is ion******************************* 
       do i=1,ninti
       zn=0
       do k=1,nepair
cc       write(3,*) i,epairlist(k,1),epairlist(k,2)
       if((i.eq.epairlist(k,1)-1000).or.(i.eq.epairlist(k,2)-1000)) zn=1
       enddo
       if(zn.eq.0) unpair=intil(i)+1000
       enddo
cc       endif
       write(*,*) 'after pairing',ninti,zn
       write(*,*) 'after unpared is ion' , unpair,ninte,nepair
!************************************************************************
!*****************determines the  time step from pair data***************
!***********************and intra atom dynamics**************************
!************************************************************************
       mindt=dt0
cc       goto 999
cc       write(4,*) 'mindt' ,mindt
       do i=1,nepair
       if((epairlist(i,2).gt.1000).and.(epairlist(i,1).lt.1000)) then
       
       dist2=ABS((xyzi(epairlist(i,2)-1000,1)-xyze(epairlist(i,1),1))**
     & 2.0+(xyzi(epairlist(i,2)-1000,2)-xyze(epairlist(i,1),2))**2.0+
     & (xyzi(epairlist(i,2)-1000,3)-xyze(epairlist(i,1),3))**2.0)
       
       tempdt=(1.0/nstep)*
     &  ( dist2**(-3.0/4.0)*ionchar(epairlist(i,2)-1000)**0.5+

     &  ( (vele(epairlist(i,1),1)-veli((epairlist(i,2)-1000),1))**2.0+
     &    (vele(epairlist(i,1),2)-veli((epairlist(i,2)-1000),2))**2.0+
     &(vele(epairlist(i,1),3)-veli((epairlist(i,2)-1000),3))**2.0 )
     &   **0.5/
     &  dist2**0.0 + omega )**(-1.0)
       
cc        write(3,*)  epairlist(i,1),epairlist(i,2)
cc       write(3,*)  tempdt,dist2

cc      write(*,*)((vele(epairlist(i,1),1)-veli((epairlist(i,2)-1000),1))
cc     & **2.0+
cc     &    (vele(epairlist(i,1),2)-veli((epairlist(i,2)-1000),2))**2.0+
cc     &(vele(epairlist(i,1),3)-veli((epairlist(i,2)-1000),3))**2.0 )**0.5,
cc      write(*,*)  'aaa'
       endif
       
       
       if((epairlist(i,2).lt.1000).and.(epairlist(i,1).lt.1000)) then
       
        dist2=ABS((xyze(epairlist(i,1),1)-xyze(epairlist(i,2),1))**2.0
     & +(xyze(epairlist(i,1),2)-xyze(epairlist(i,2),2))**2.0+
     & (xyze(epairlist(i,1),3)-xyze(epairlist(i,2),3))**2.0)**(1.0/2.0)
       
cc        write(3,*) dist2
        tempdt=(1.0/nstep)*
     &  (dist2**(-3.0/4.0)+
     &  ((vele(epairlist(i,1),1)-vele(epairlist(i,2),1))**2.0+
     &   (vele(epairlist(i,1),2)-vele(epairlist(i,2),2))**2.0+
     &   (vele(epairlist(i,1),3)-vele(epairlist(i,2),3))**2.0)**0.5/
     &  dist2**0.5*0.0 + omega)**(-1.0)
cc       write(3,*)  epairlist(i,1),epairlist(i,2)         
cc       write(3,*)  tempdt,dist2
cc       write(*,*) '******'
cc       write(*,*) xyze(epairlist(i,1),1), xyze(epairlist(i,2),1)
cc       write(*,*) xyze(epairlist(i,1),2), xyze(epairlist(i,2),2)
cc       write(*,*) xyze(epairlist(i,1),3), xyze(epairlist(i,2),3) 
       endif       
       if((epairlist(i,2).lt.1000).and.(epairlist(i,1).gt.1000)) then
       
        dist2=ABS((xyzi(epairlist(i,1)-1000,1)-xyze(epairlist(i,2),1))**2.0
     &+(xyzi(epairlist(i,1)-1000,2)-xyze(epairlist(i,2),2))**2.0+
     &  (xyzi(epairlist(i,1)-1000,3)-xyze(epairlist(i,2),3))**2.0)
       

        tempdt=(1.0/nstep)*
     &  (dist2**(-3.0/4.0)*ionchar(epairlist(i,1)-1000)**0.5+
     &  ((vele(epairlist(i,2),1)-veli((epairlist(i,1)-1000),1))**2.0+
     &  (vele(epairlist(i,2),2)-veli((epairlist(i,1)-1000),2))**2.0+
     &  (vele(epairlist(i,2),3)-veli((epairlist(i,1)-1000),3))**2.0)
     &   **0.5/
     &  dist2**0.0 + omega)**(-1.0)
 
      write(*,*)((vele(epairlist(i,1),1)-veli((epairlist(i,2)-1000),1))
     &**2.0+
     &    (vele(epairlist(i,1),2)-veli((epairlist(i,2)-1000),2))**2.0+
     &(vele(epairlist(i,1),3)-veli((epairlist(i,2)-1000),3))**2.0 )**0.5
       write(*,*) 'aaa'
       endif
       
       if((epairlist(i,2).gt.1000).and.(epairlist(i,1).gt.1000)) then
       
       
       dist2=ABS((xyzi(epairlist(i,1)-1000,1)-
     &  xyzi(epairlist(i,2)-1000,1))**2.0+
     &  (xyzi(epairlist(i,1)-1000,2)-xyzi(epairlist(i,2)-1000,2))**2.0+
     &  (xyzi(epairlist(i,1)-1000,3)-xyzi(epairlist(i,2)-1000,3))**2.0)
     & **(1.0/2.0)
    
       tempdt=(1.0/nstep)*
     &  ( dist2**(-3.0/4.0)*ionchar(epairlist(i,2)-1000)**0.5*
     &  ionchar(epairlist(i,1)-1000)**0.5/mi**0.5+

     & ((veli(epairlist(i,1)-1000,1)-veli(epairlist(i,2)-1000,1))**2.0+
     &  (veli(epairlist(i,1)-1000,2)-veli(epairlist(i,2)-1000,2))**2.0+
     &  (veli(epairlist(i,1)-1000,3)-veli(epairlist(i,2)-1000,3))**2.0 )
     &  **0.5/
     &  dist2**0.5 + omega)**(-1.0)
              
       endif
       if (tempdt.lt.mindt) mindt=tempdt
cc       write(3,*) mindt,epairlist(i,1),epairlist(i,2),dist2,i,'**'
       enddo
cc       write(*,*) 'before assigned electrons', tempdt

!**************************checking assigned electrons***************************
!          do i=1,512
           do i=1,nele
	   if(ionhist(i).eq.0) then 
	   dist2=ABS((xyzi(eiasm(i),1)-xyze(i,1))**
     & 2.0+(xyzi(eiasm(i),2)-xyze(i,2))**2.0+
     & (xyzi(eiasm(i),3)-xyze(i,3))**2.0)
       
       if(dist2.gt.0.0) then
       tempdt=(1.0/nstep)*
     &  (dist2**(-3.0/4.0)*ifeiint0(i,eiasm(i))**0.5+
     &  ((vele(i,1)-veli(eiasm(i),1))**2.0+
     &   (vele(i,2)-veli(eiasm(i),2))**2.0+
     &   (vele(i,3)-veli(eiasm(i),3))**2.0)**0.5/
     &  dist2**0.5 + omega)**(-1.0)
	   endif
	   endif 
	   if (tempdt.lt.mindt) mindt=tempdt
cc           write(*,*) tempdt, 'tempdt',dist2**0.5
cc        write(3,*) mindt
	   enddo
cc        write(*,*) 'after assigned electrons', tempdt,mindt
	   
       write(*,*) 'before log' 
       write(*,*) log(dt0/mindt),mindt
       nmax=(1.0+Int(log(dt0/mindt)/log(2.0)))
       if(nmax.gt.6.0) nmax=6.0
cc       nmax=12.0
cc       write(*,*) nmax
       tempdt=dt0/2.0**nmax
cc        write(*,*) tempdt, 'tempdt'
cc       write(4,*) 'tempdt', tempdt
      write(*,*) 'nmax' , nmax
!****************************sets grid in a time domain**************************
       if(tempdt.le.dt) then 
       dt=tempdt
cc       write(*,*) 'dt=tempdt'
       else
       
       if (ABS((tempt/tempdt)-nINT(tempt/tempdt)).lt.1.0/2.0**64) 
     &  dt=tempdt 
cc        dt=tempdt
        write(*,*) 'changing up', tempt/tempdt,tempdt,dt
        write(*,*) 'after'
        write(*,*) ABS((tempt/tempdt)-nINT(tempt/tempdt))
        write(*,*) 1.0/2.0**64
        write(*,*) 'after 1'
cc        endif

       endif
cc999    dt=dt0       
       write(*,*) dt0/dt
       tempt=tempt+dt
       if(tempt.gt.dt0) tempt=0.0
 
cc       write(*,*) 'temporary time', tempt,dt
 
!*********************************big endif here*********************************
       endif
!*************Solves Coulomb problem for bound  electron  with host ion**********
!***Solves Coulomb problem for  ionized      electron  with strongest  ion*******
!*******************or strongest electron-electron and the laser*****************
!**********original positions and velocities for parametric integration**********
       do i=1,nion
       xyzio(i,1)=xyzi(i,1)
       xyzio(i,2)=xyzi(i,2)
       xyzio(i,3)=xyzi(i,3)
       velio(i,1)=veli(i,1)
       velio(i,2)=veli(i,2)
       velio(i,3)=veli(i,3)
       enddo
!*******************ionized with ion or e-e and laser****************************
!********************or unpaired with laser if one left**************************
       do i=1,nepair
       if((epairlist(i,2).gt.1000).and.(epairlist(i,1).lt.1000)) then
       m1=me
       m2=mi
	   kk=ionchar(epairlist(i,2)-1000)
cc       write(*,*) kk, 'e-i *************************************'
       mu=m1*m2/(m1+m2)
       xx=-xyze(epairlist(i,1),1)+xyzi(epairlist(i,2)-1000,1)
       yy=-xyze(epairlist(i,1),2)+xyzi(epairlist(i,2)-1000,2)
       zz=-xyze(epairlist(i,1),3)+xyzi(epairlist(i,2)-1000,3)
       px=-vele(epairlist(i,1),1)+veli(epairlist(i,2)-1000,1)
       py=-vele(epairlist(i,1),2)+veli(epairlist(i,2)-1000,2)
       pz=-vele(epairlist(i,1),3)+veli(epairlist(i,2)-1000,3)
       Rxx=(m1*xyze(epairlist(i,1),1)+
     & m2*xyzi(epairlist(i,2)-1000,1))/(m1+m2)
       Ryy=(m1*xyze(epairlist(i,1),2)+
     & m2*xyzi(epairlist(i,2)-1000,2))/(m1+m2)
       Rzz=(m1*xyze(epairlist(i,1),3)+
     & m2*xyzi(epairlist(i,2)-1000,3))/(m1+m2)
       Vxx=(m1*vele(epairlist(i,1),1)+
     & m2*veli(epairlist(i,2)-1000,1))/(m1+m2)
       Vyy=(m1*vele(epairlist(i,1),2)+
     & m2*veli(epairlist(i,2)-1000,2))/(m1+m2)
       Vzz=(m1*vele(epairlist(i,1),3)+
     & m2*veli(epairlist(i,2)-1000,3))/(m1+m2)
cc       write(3,*) 'before calling coulomb',dt , 'dt just before'
cc       write(3,*) xx,'xx'
cc       write(3,*) yy,'yy'
cc       write(3,*) zz,'zz'
cc       write(3,*) px,'px'
cc       write(3,*) py,'py'
cc       write(3,*) pz,'pz'
cc       write(3,*) Rxx,'Rxx'
cc       write(3,*) Ryy,'Ryy'
cc       write(3,*) Rzz,'Rzz'
cc       write(3,*) Vxx,'Vxx'
cc       write(3,*) Vyy,'Vyy'
cc       write(3,*) Vzz,'Vzz'
       Call Coulomb(xx,yy,zz,px,py,pz,dt,mu,kk)
       Call Volcov(Rxx,Ryy,Rzz,Vxx,Vyy,Vzz,tt,dt,ww,mu,E00,kk)
cc       write(3,*) 'after calling coulomb-volcov',xx,yy,zz,px,py,pz
       xyze(epairlist(i,1),1)=Rxx-m2/(m1+m2)*xx
       xyze(epairlist(i,1),2)=Ryy-m2/(m1+m2)*yy
       xyze(epairlist(i,1),3)=Rzz-m2/(m1+m2)*zz
       xyzi(epairlist(i,2)-1000,1)=Rxx+m1/(m1+m2)*xx
       xyzi(epairlist(i,2)-1000,2)=Ryy+m1/(m1+m2)*yy
       xyzi(epairlist(i,2)-1000,3)=Rzz+m1/(m1+m2)*zz
       vele(epairlist(i,1),1)=Vxx-m2/(m1+m2)*px
       vele(epairlist(i,1),2)=Vyy-m2/(m1+m2)*py
       vele(epairlist(i,1),3)=Vzz-m2/(m1+m2)*pz
       veli(epairlist(i,2)-1000,1)=Vxx+m1/(m1+m2)*px
       veli(epairlist(i,2)-1000,2)=Vyy+m1/(m1+m2)*py
       veli(epairlist(i,2)-1000,3)=Vzz+m1/(m1+m2)*pz
       endif
       
       if((epairlist(i,2).lt.1000).and.(epairlist(i,1).lt.1000)) then
       m1=me
       m2=me
	   kk=-1.0
cc       write(*,*) kk, 'e-e *************************************'
       mu=m1*m2/(m1+m2)
       xx=-xyze(epairlist(i,1),1)+xyze(epairlist(i,2),1)
       yy=-xyze(epairlist(i,1),2)+xyze(epairlist(i,2),2)
       zz=-xyze(epairlist(i,1),3)+xyze(epairlist(i,2),3)
       px=-vele(epairlist(i,1),1)+vele(epairlist(i,2),1)
       py=-vele(epairlist(i,1),2)+vele(epairlist(i,2),2)
       pz=-vele(epairlist(i,1),3)+vele(epairlist(i,2),3)
       Rxx=(m1*xyze(epairlist(i,1),1)+m2*xyze(epairlist(i,2),1))/(m1+m2)
       Ryy=(m1*xyze(epairlist(i,1),2)+m2*xyze(epairlist(i,2),2))/(m1+m2)
       Rzz=(m1*xyze(epairlist(i,1),3)+m2*xyze(epairlist(i,2),3))/(m1+m2)
       Vxx=(m1*vele(epairlist(i,1),1)+m2*vele(epairlist(i,2),1))/(m1+m2)
       Vyy=(m1*vele(epairlist(i,1),2)+m2*vele(epairlist(i,2),2))/(m1+m2)
       Vzz=(m1*vele(epairlist(i,1),3)+m2*vele(epairlist(i,2),3))/(m1+m2)
cc       write(3,*) 'before calling coulomb',dt , 'dt just before'
cc       write(3,*) xx,'xx'
cc       write(3,*) yy,'yy'
cc       write(3,*) zz,'zz'
cc       write(3,*) px,'px'
cc       write(3,*) py,'py'
cc       write(3,*) pz,'pz'
cc       write(3,*) Rxx,'Rxx'
cc       write(3,*) Ryy,'Ryy'
cc       write(3,*) Rzz,'Rzz'
cc       write(3,*) Vxx,'Vxx'
cc       write(3,*) Vyy,'Vyy'
cc       write(3,*) Vzz,'Vzz'
       Call Coulomb(xx,yy,zz,px,py,pz,dt,mu,kk)
       Call Volcov(Rxx,Ryy,Rzz,Vxx,Vyy,Vzz,tt,dt,ww,2.0*me,E0,-kk)
cc       write(3,*) 'after calling coulomb-volcov',xx,yy,zz,px,py,pz
       xyze(epairlist(i,1),1)=Rxx-m2/(m1+m2)*xx
       xyze(epairlist(i,1),2)=Ryy-m2/(m1+m2)*yy
       xyze(epairlist(i,1),3)=Rzz-m2/(m1+m2)*zz
       xyze(epairlist(i,2),1)=Rxx+m1/(m1+m2)*xx
       xyze(epairlist(i,2),2)=Ryy+m1/(m1+m2)*yy
       xyze(epairlist(i,2),3)=Rzz+m1/(m1+m2)*zz
       vele(epairlist(i,1),1)=Vxx-m2/(m1+m2)*px
       vele(epairlist(i,1),2)=Vyy-m2/(m1+m2)*py
       vele(epairlist(i,1),3)=Vzz-m2/(m1+m2)*pz
       vele(epairlist(i,2),1)=Vxx+m1/(m1+m2)*px
       vele(epairlist(i,2),2)=Vyy+m1/(m1+m2)*py
       vele(epairlist(i,2),3)=Vzz+m1/(m1+m2)*pz
       
       endif
       
       if((epairlist(i,2).lt.1000).and.(epairlist(i,1).gt.1000)) then
       m1=mi
       m2=me
       mu=m1*m2/(m1+m2)
       kk=ionchar(epairlist(i,1)-1000)
cc       write(*,*) kk, 'i-e *************************************'
       xx=-xyzi(epairlist(i,1)-1000,1)+xyze(epairlist(i,2),1)
       yy=-xyzi(epairlist(i,1)-1000,2)+xyze(epairlist(i,2),2)
       zz=-xyzi(epairlist(i,1)-1000,3)+xyze(epairlist(i,2),3)
       px=-veli(epairlist(i,1)-1000,1)+vele(epairlist(i,2),1)
       py=-veli(epairlist(i,1)-1000,2)+vele(epairlist(i,2),2)
       pz=-veli(epairlist(i,1)-1000,3)+vele(epairlist(i,2),3)
       Rxx=(m1*xyzi(epairlist(i,1)-1000,1)
     & +m2*xyze(epairlist(i,2),1))/(m1+m2)
       Ryy=(m1*xyzi(epairlist(i,1)-1000,2)
     & +m2*xyze(epairlist(i,2),2))/(m1+m2)
       Rzz=(m1*xyzi(epairlist(i,1)-1000,3)+
     & m2*xyze(epairlist(i,2),3))/(m1+m2)
       Vxx=(m1*veli(epairlist(i,1)-1000,1)+
     & m2*vele(epairlist(i,2),1))/(m1+m2)
       Vyy=(m1*veli(epairlist(i,1)-1000,2)+
     & m2*vele(epairlist(i,2),2))/(m1+m2)
       Vzz=(m1*veli(epairlist(i,1)-1000,3)+
     & m2*vele(epairlist(i,2),3))/(m1+m2)
cc       write(3,*) 'before calling coulomb',dt , 'dt just before'
cc       write(3,*) xx,'xx'
cc       write(3,*) yy,'yy'
cc       write(3,*) zz,'zz'
cc       write(3,*) px,'px'
cc       write(3,*) py,'py'
cc       write(3,*) pz,'pz'
cc       write(3,*) Rxx,'Rxx'
cc       write(3,*) Ryy,'Ryy'
cc       write(3,*) Rzz,'Rzz'
cc       write(3,*) Vxx,'Vxx'
cc       write(3,*) Vyy,'Vyy'
cc       write(3,*) Vzz,'Vzz'
       Call Coulomb(xx,yy,zz,px,py,pz,dt,mu,kk)
       Call Volcov(Rxx,Ryy,Rzz,Vxx,Vyy,Vzz,tt,dt,ww,mu,E00,kk)
cc       write(3,*) 'after calling coulomb-volcov',xx,yy,zz,px,py,pz
       xyzi(epairlist(i,1)-1000,1)=Rxx-m2/(m1+m2)*xx
       xyzi(epairlist(i,1)-1000,2)=Ryy-m2/(m1+m2)*yy
       xyzi(epairlist(i,1)-1000,3)=Rzz-m2/(m1+m2)*zz
       xyze(epairlist(i,2),1)=Rxx+m1/(m1+m2)*xx
       xyze(epairlist(i,2),2)=Ryy+m1/(m1+m2)*yy
       xyze(epairlist(i,2),3)=Rzz+m1/(m1+m2)*zz
       veli(epairlist(i,1)-1000,1)=Vxx-m2/(m1+m2)*px
       veli(epairlist(i,1)-1000,2)=Vyy-m2/(m1+m2)*py
       veli(epairlist(i,1)-1000,3)=Vzz-m2/(m1+m2)*pz
       vele(epairlist(i,2),1)=Vxx+m1/(m1+m2)*px
       vele(epairlist(i,2),2)=Vyy+m1/(m1+m2)*py
       vele(epairlist(i,2),3)=Vzz+m1/(m1+m2)*pz
       endif
       if ((epairlist(i,2).gt.1000).and.(epairlist(i,1).gt.1000)) then
       m1=mi
       m2=mi
       mu=m1*m2/(m1+m2)
       kk=-ionchar(epairlist(i,1)-1000)*ionchar(epairlist(i,2)-1000)
cc       write(*,*) kk, 'i-i *************************************'
cc       write(*,*) 'aaaa'
cc       write(*,*) epairlist(i,1), epairlist(i,2)
cc       write(*,*) ionchar(epairlist(i,1)-1000),
cc     &                      ionchar(epairlist(i,2)-1000)
cc       write(*,*) 'bbbb'
       xx=-xyzi(epairlist(i,1)-1000,1)+xyzi(epairlist(i,2)-1000,1)
       yy=-xyzi(epairlist(i,1)-1000,2)+xyzi(epairlist(i,2)-1000,2)
       zz=-xyzi(epairlist(i,1)-1000,3)+xyzi(epairlist(i,2)-1000,3)
       px=-veli(epairlist(i,1)-1000,1)+veli(epairlist(i,2)-1000,1)
       py=-veli(epairlist(i,1)-1000,2)+veli(epairlist(i,2)-1000,2)
       pz=-veli(epairlist(i,1)-1000,3)+veli(epairlist(i,2)-1000,3)
       Rxx=(m1*xyzi(epairlist(i,1)-1000,1)+
     & m2*xyzi(epairlist(i,2)-1000,1))/(m1+m2)
       Ryy=(m1*xyzi(epairlist(i,1)-1000,2)+
     & m2*xyzi(epairlist(i,2)-1000,2))/(m1+m2)
       Rzz=(m1*xyzi(epairlist(i,1)-1000,3)+
     & m2*xyzi(epairlist(i,2)-1000,3))/(m1+m2)
       Vxx=(m1*veli(epairlist(i,1)-1000,1)+
     & m2*veli(epairlist(i,2)-1000,1))/(m1+m2)
       Vyy=(m1*veli(epairlist(i,1)-1000,2)+
     & m2*veli(epairlist(i,2)-1000,2))/(m1+m2)
       Vzz=(m1*veli(epairlist(i,1)-1000,3)+
     & m2*veli(epairlist(i,2)-1000,3))/(m1+m2)
cc       write(3,*) xx,yy,zz,px,py,pz,'before c ii'
cc       write(3,*) epairlist(i,1),epairlist(i,2)
       Call Coulomb(xx,yy,zz,px,py,pz,dt,mu,kk)
cc       write(3,*) xx,yy,zz,px,py,pz,'after  c ii'
cc       stop
!*********Center free evolution***************************************
       Call Volcov(Rxx,Ryy,Rzz,Vxx,Vyy,Vzz,tt,dt,ww,mu,E00,kk)
       xyzi(epairlist(i,1)-1000,1)=Rxx-m2/(m1+m2)*xx
       xyzi(epairlist(i,1)-1000,2)=Ryy-m2/(m1+m2)*yy
       xyzi(epairlist(i,1)-1000,3)=Rzz-m2/(m1+m2)*zz
       xyzi(epairlist(i,2)-1000,1)=Rxx+m1/(m1+m2)*xx
       xyzi(epairlist(i,2)-1000,2)=Ryy+m1/(m1+m2)*yy
       xyzi(epairlist(i,2)-1000,3)=Rzz+m1/(m1+m2)*zz
       veli(epairlist(i,1)-1000,1)=Vxx-m2/(m1+m2)*px
       veli(epairlist(i,1)-1000,2)=Vyy-m2/(m1+m2)*py
       veli(epairlist(i,1)-1000,3)=Vzz-m2/(m1+m2)*pz
       veli(epairlist(i,2)-1000,1)=Vxx+m1/(m1+m2)*px
       veli(epairlist(i,2)-1000,2)=Vyy+m1/(m1+m2)*py
       veli(epairlist(i,2)-1000,3)=Vzz+m1/(m1+m2)*pz
       endif
       
       
       enddo
cc       stop
cc       write(*,*) 'after paired',unpair
!****************Treats a single unpaired particle if any************************  
!***************************as free in a laser***********************************   
cc       unpair=0 
       if (unpair.gt.0)  then
	   i=unpair
cc       write(*,*) 'am i here'
       if (unpair.gt.1000) then      
       mu=mi
       xx=xyzi(i-1000,1)
       yy=xyzi(i-1000,2)
       zz=xyzi(i-1000,3)
       px=veli(i-1000,1)
       py=veli(i-1000,2)
       pz=veli(i-1000,3)
       kk=-corechar
       Call Volcov(xx,yy,zz,px,py,pz,tt,dt,ww,mu,E0,kk)
       xyzi(i-1000,1)=xx
       xyzi(i-1000,2)=yy
       xyzi(i-1000,3)=zz
       veli(i-1000,1)=px
       veli(i-1000,2)=py
       veli(i-1000,3)=pz
       else  
cc       write(*,*) 'why ?' ,unpair
       mu=me
       xx=xyze(i,1)
       yy=xyze(i,2)
       zz=xyze(i,3)
       px=vele(i,1)
       py=vele(i,2)
       pz=vele(i,3)  
       kk=1.0
cc       write(*,*) 'before volcov',xx,yy,zz 
       Call Volcov(xx,yy,zz,px,py,pz,tt,dt,ww,mu,E0,kk)
       xyze(i,1)=xx
       xyze(i,2)=yy
       xyze(i,3)=zz
       vele(i,1)=px
       vele(i,2)=py
       vele(i,3)=pz
cc       write(*,*) 'what am i doing here'
       endif      
       endif
cc       write(*,*) 'before return'
cc       return 

cc       write(3,*) 'next after integrating'
cc       do i=1,ninti
cc       write(3,*) xyzi(intil(i),1),xyzi(intil(i),2),xyzi(intil(i),3)
cc       enddo
cc       write(3,*) 'elec1'
cc       do i=1,ninte
cc       write(3,*) xyze(intel(i),1),xyze(intel(i),2),xyze(intel(i),3)
cc       enddo
      

!************************bound with host ****************************************
!**********this has been moved because of parametric integration*****************
!       do i=1,512
       do i=1,nele
       if((ionhist(i).eq.0).or.(ionhist(i).eq.2)) then
       xx=xyze(i,1)-xyzio(eiasm(i),1)
       yy=xyze(i,2)-xyzio(eiasm(i),2)
       zz=xyze(i,3)-xyzio(eiasm(i),3)
       px=vele(i,1)-velio(eiasm(i),1)
       py=vele(i,2)-velio(eiasm(i),2)
       pz=vele(i,3)-velio(eiasm(i),3)
       mu=me
       kk=ifeiint0(i,eiasm(i))
cc       kk=1.0
cc       write(*,*) 'before calling coulomb',unpair,dt,'dt'
cc       write(*,*) kk, 'bwh  *************************************'
cc       write(*,*) xx,yy,zz
cc       write(*,*) px,py,pz
cc       write(*,*) 'before'
cc       write(*,*) dt
       Call Coulomb(xx,yy,zz,px,py,pz,dt,mu,kk)
cc       write(*,*) 'after'
cc       write(*,*) xx,yy,zz
cc       write(*,*) px,py,pz
       vele(i,1)=px+veli(eiasm(i),1)
       vele(i,2)=py+veli(eiasm(i),2)
       vele(i,3)=pz+veli(eiasm(i),3)
       xyze(i,1)=xx+xyzi(eiasm(i),1)
       xyze(i,2)=yy+xyzi(eiasm(i),2)
       xyze(i,3)=zz+xyzi(eiasm(i),3)
       endif
       enddo
       do i=1,nion
       if (ionchar(i).lt.1.0) then
       xyzi(i,1)=xyzi(i,1)+veli(i,1)*dt
       xyzi(i,2)=xyzi(i,2)+veli(i,2)*dt
       xyzi(i,3)=xyzi(i,3)+veli(i,3)*dt
       endif
       enddo
cc       write(*,*) 'after assigned',unpair
 
       return 
       end

       
       Subroutine Coulomb(xx,yy,zz,px,py,pz,dt,mu,kk)
	   real*8 mu
	   real*8 kk
	   real*8 llx
	   real*8 lly
	   real*8  llz
	   real*8 px
	   real*8 py
	   real*8 pz
	   real*8 xx
	   real*8 yy
	   real*8 zz
	   real*8 rrr
	   real*8 lll
	   real*8 elx
	   real*8 ely
	   real*8 elz
	   real*8 erx
	   real*8 ery
	   real*8 erz
	   real*8 e3x
	   real*8 e3y
	   real*8 e3z
	   real*8 xxx
	   real*8 yyy
	   real*8 zzz
	   real*8 pxx
	   real*8 pyy
	   real*8 pzz
	   real*8 EE
	   real*8 a
	   real*8 tt0
	   real*8 tt1
	   real*8 int1
	   real*8 x1
	   real*8 x2
	   real*8 dtheta
	   real*8 fmid
	   real*8 fff
	   real*8 tt
	   real*8 dx
	   real*8 xmid
	   real*8 ttacc
	   real*8 rtt
	   real*8 xxxxtt
	   real*8 yyyytt
	   real*8 xxxxxtt
	   real*8 yyyyytt
	   real*8 zzzzztt
	   real*8 veltt
	   real*8 velrr
	   real*8 errx
	   real*8 erry
	   real*8 ettx
	   real*8 etty
	   real*8 vxxxxtt
	   real*8 vyyyytt
	   real*8 vxxxxxtt
	   real*8 vyyyyytt
	   real*8 vzzzzztt
	   real*8 one
           real*8 aracos
           real*8 aa
           real*8 r0 
           real*8 vv0 
           real*8 svv0 
           real*8 asv
           real*8 omegak
           real*8 vv 
           real*8 ttt
           real*8 dt
           real*8 vv1
           real*8 d1,d2
           integer zn
cc           write(*,*) 'i am in coulomb',dt
	   dtheta=3.141/2.0
	   ttacc=0.001*dt 
           one=1.0
	   if(kk.lt.0.0) one=-1.0
c           px=0.0
c          py=0.4
c           pz=0.0
c           xx=1.0
c           yy=0.0
c           zz=0.0 
	   
       llx= ( pz*yy-py*zz)*mu
       lly= (-pz*xx+px*zz)*mu
       llz= ( py*xx-px*yy)*mu
c       write(*,*) llx,'llx'
c       write(*,*) lly,'lly'
c       write(*,*) llz,'llz'
       rrr=(xx**2.0+yy**2.0+zz**2.0)**0.5
       lll=(llx**2.0+lly**2.0+llz**2.0)**0.5
cc       write(7,*) 'before',xx,yy,zz,px,py,pz,'before'
!**********plain leap frog sub-jump for zero angular momentum************************
        if((lll.eq.0.0).or.(mu.gt.500.0).or.(rrr.gt.3000.0)) then 
cc        if(lll.ne.1.11111) then
cc         if(rrr.lt.0.1) then
cc       if (mu.gt.200) then
cc       write(7,*) 'before',xx,yy,zz,px,py,pz,'before'
cc       endif
 
cc       write(*,*) 'am i in leap ?'
cc       write(1,*) lll 
       xx=xx+0.5*dt*px
       yy=yy+0.5*dt*py
       zz=zz+0.5*dt*pz
       rrr=(xx**2.0+yy**2.0+zz**2.0+0.000001)**0.5
       px=px-dt/mu*xx*kk/rrr**(3.0)
       py=py-dt/mu*yy*kk/rrr**(3.0)
       pz=pz-dt/mu*zz*kk/rrr**(3.0)
       xx=xx+0.5*dt*px
       yy=yy+0.5*dt*py
       zz=zz+0.5*dt*pz
cc       if (mu.gt.200) then 
cc       write(3,*) 'after',xx,yy,zz,px,py,pz
cc       endif
       return
       endif
       elx=llx/lll
       ely=lly/lll
       elz=llz/lll
       erx=xx/rrr
       ery=yy/rrr
       erz=zz/rrr
       e3x= elz*ery-ely*erz
       e3y=-elz*erx+elx*erz
       e3z= ely*erx-elx*ery
       xxx=erx*xx+ery*yy+erz*zz
       yyy=e3x*xx+e3y*yy+e3z*zz
       zzz=elx*xx+ely*yy+elz*zz
       pxx=erx*px+ery*py+erz*pz
       pyy=e3x*px+e3y*py+e3z*pz
       pzz=elx*px+ely*py+elz*pz
       lll=(pyy*xxx-pxx*yyy)*mu
       EE=mu*pxx**2.0/2.0+mu*pyy**2.0/2.0 -kk/rrr
       a=(Abs(1.0+ 2.0*EE*lll**2/(mu*kk**2.0)))**0.5
cc       write(7,*) a,EE,lll,'aEElll'

cc       if((ABS(a-1.0).lt.0.00000000000000001).and.
cc     & (EE.ne.0.0)) then
       if((a.eq.1.0).and.(EE.ne.0.0)) then
       write(*,*) 'am i in leap 2?'
       xx=xx+0.5*dt*px
       yy=yy+0.5*dt*py
       zz=zz+0.5*dt*pz
       rrr=(xx**2.0+yy**2.0+zz**2.0)**0.5
       px=px-dt/mu*xx*kk/rrr**(3.0)
       py=py-dt/mu*yy*kk/rrr**(3.0)
       pz=pz-dt/mu*zz*kk/rrr**(3.0)
       xx=xx+0.5*dt*px
       yy=yy+0.5*dt*py
       zz=zz+0.5*dt*pz
       return
       endif

cc       write(*,*) a,mu,lll,kk,EE, 'a before'
!********note this ArcTan is opposite to Mathematica !!!************************
       tt0=Atan2(yyy,xxx)
       if(a.gt.0.0) then
       aracos=(lll**2.0/(rrr*mu*ABS(kk))-one)/a
cc       write(*,*) 'aracos before Kepler', aracos,ACos(aracos) 
       if(aracos.ge.1.0) aracos=1.0
       if(aracos.le.-1.0) aracos=-1.0
c       tt1=tt0-ACos((lll**2.0/(rrr*mu*ABS(kk))-one)/a)
c     & *sign(1.0,pxx*pyy) 
        tt1=tt0-ACos(aracos)
     & *sign(1.0,pxx*pyy)
cc        write(*,*) 'angles before Kepler',tt0,tt1,aracos,ACos(aracos)
!*****************solves Kepler angle vv versus dt equation**********************
!**************************using iteration method********************************
!************************this should be well tested******************************
!************************the case of parabola ***********************************
c       dt=3.141*4.0/1.0/100.0
c       dt=0.0
       if((a.eq.1.0).and.(EE.eq.0)) then 
       aa=tan(-tt1/2.0)
       aa=dt*mu*kk**2.0/lll**2.0 +aa+aa**3.0
       tt=(54.0*aa+(864.0+2916.0*aa**2.0))**(1.0/3.0)
       tt=2.0*2.0**(1.0/3.0)/tt+tt/(3.0*2.0**(1.0/3.0))
       tt=tt1+2.0*atan(tt)
       rtt=(lll**2.0/mu/ABS(kk))/(one+a*Cos(tt-tt1))
cc       write(*,*) 'am i in parabola ?'
       endif
!**************************the case of hiperbola*********************************
       if((a.gt.1.0).and.(EE.lt.0.0)) write(1,*) aa,EE
       if((a.lt.1.0).and.(EE.gt.0.0)) write(1,*) aa,EE
       if(a.gt.1.0) then 
cc       write(*,*) a, 'i am in hiperbola',mu
       r0=kk/(2.0*EE)
       r0=ABS(r0)
       omegak=(ABS(kk)/mu/r0**3.0)**0.5*sign(1.0,lll)
cc       write(*,*) omegak*dt,a
cc       write(*,*) omegak,r0,EE,sign(1.0,lll),'aaaa'
       vv0=(r0*one+rrr)/a/r0
       if(vv0.lt.1.0) vv0=1.0
cc       write(*,*) vv0,'vv0',r0
       vv0=log(vv0+(vv0**2.0-1.0)**0.5)*sign(1.0,pxx*pyy)
       svv0=sinh(vv0)
       vv=vv0
       zn=0
10     continue
cc       do i=1,200
       vv1=vv
       asv=dt*omegak/a+svv0+(vv-vv0)/a*one
       vv=log(asv+(1.0+asv**2.0)**0.5)
cc       write(2,*) vv,'vv'
cc       enddo
cc       if(zn.eq.1) d1=ABS(vv1-vv)
cc       if(zn.eq.-1) d2=ABS(vv1-vv)
cc       if(ABS(dd1-dd2).lt.000000000000001) goto 11
cc       zn=-zn
       zn=zn+1 
       if((ABS(vv1-vv).le.0.00000001).or.(zn.gt.10000)) goto 11
       goto 10
11     continue
       rtt=r0*(a*cosh(vv)-one)
       endif
!**************************the case of ellipse***********************************
       if(a.lt.1.0) then
cc       if(a.gt.0.99) a=0.99
       r0=-kk/(2.0*EE)
cc       write(7,*) 'i am in ellipse',a
cc       write(*,*) r0
cc       write(*,*) kk
cc       write(*,*) EE
       omegak=(ABS(kk)/mu/r0**3.0)**0.5*sign(1.0,lll)
cc       write(*,*) omegak*dt,a
cc       write(*,*) 2*3.14/omegak,'omegac'
       aracos=(r0-rrr)/a/r0
cc       write(7,*) 'warning aracos !!!', aracos
       if(aracos.ge.1.0) aracos=1.00
       if(aracos.le.-1.0) aracos=-1.00
       vv0=Acos(aracos)*sign(1.0,pxx*pyy)
       svv0=a*sin(vv0)
       vv=vv0 
       zn=0
20     continue
cc       do i=1,200
       vv1=vv
       vv=vv0+dt*omegak-svv0+a*sin(vv)
cc       write(7,*) tt1, 'vv',2*3.14/omegak,dt
cc       write(*,*) vv
cc       enddo
cc       if(zn.eq.1) d1=ABS(vv1-vv)
cc       if(zn.eq.-1) d2=ABS(vv1-vv)
cc       if(ABS(d1-dd2).lt.0.000000000000001) goto 21
cc       zn=-zn
       zn=zn+1
       if((ABS(vv1-vv).le.0.00000001).or.(zn.gt.10000)) goto 21 
       goto 20
21     continue
       rtt=r0*(1.0-a*cos(vv))
       endif
cc       write(7,*) 'after ellipse', vv,rtt
!********************************************************************************
       if(a.ne.0.0) then 
       aracos=(lll**2.0/(rtt*mu*ABS(kk))-one)/a
       if(aracos.gt.1.0) aracos=1.0
       if(aracos.lt.-1.0) aracos=-1.0
       tt=tt1+ACos(aracos)*sign(1.0,pxx*pyy)
       if (sign(1.0,tt)*sign(1.0,lll).lt.0) 
     & tt=tt1-ACos(aracos)*sign(1.0,pxx*pyy)
       endif 
c       tt=t0+dt
cc       write(7,*) 'angle ', tt,tt1,aracos
!********************************************************************************* 
       else
c       dt=3.141*2.0/1.0/100.0
       tt=(mu*kk**2.0)/lll**3.0*dt+tt0
cc       tt=dt
       tt1=tt
cc       dt=3.141*2.0/2.0/10.0*sign(1.0,lll)
cc       tt=dt+tt0
       endif 
c       write(1,*) tt,dt
c       dt=3.141*2.0/2.0/100.0*sign(1.0,lll)
c       tt=dt
       rtt=(lll**2.0/mu/ABS(kk))/(one+a*Cos(tt-tt1))  
       xxxxtt=rtt*Cos(tt)
       yyyytt=rtt*Sin(tt)
       xxxxxtt=(erx*xxxxtt+e3x*yyyytt)
       yyyyytt=(ery*xxxxtt+e3y*yyyytt)
       zzzzztt=(erz*xxxxtt+e3z*yyyytt)
       xx=xxxxxtt
       yy=yyyyytt
       zz=zzzzztt
       veltt=lll/rtt/mu
       velrr=(2.0/mu)**0.5*
     & abs(kk/rtt-lll**2.0/(2.0*rtt**2.0*mu)+EE)**0.5
       velrr=velrr*sign(1.0,Sin(tt-tt1)*lll)
       errx= Cos(tt)
       erry= Sin(tt)
       ettx=-Sin(tt)
       etty= Cos(tt)
       vxxxxtt=velrr*errx + veltt*ettx
       vyyyytt=velrr*erry + veltt*etty
       vxxxxxtt=(erx*vxxxxtt+e3x*vyyyytt)
       vyyyyytt=(ery*vxxxxtt+e3y*vyyyytt)
       vzzzzztt=(erz*vxxxxtt+e3z*vyyyytt)
       px=vxxxxxtt
       py=vyyyyytt
       pz=vzzzzztt
cc       write(*,*) mu*veltt**2.0/2.0+mu*velrr**2.0/2.0 -kk/rtt 
cc       write(*,*) 'is conserved ?', EE
cc       write(*,*) a
cc       write(7,*) 'leaving coulomb',xx,yy,zz,px,py,pz
cc       write(*,*) 'leaving coulomb',tt
cc       write(*,*) 'leaving coulomb', mu, rtt,kk,lll,EE,velrr
cc       write(*,*) 
cc     & abs( (kk/rtt-lll**2.0/(2.0*rtt**2.0*mu)+EE))**0.5
cc       write(*,*) 'bye'
cc       write(*,*) 'lc',xx
cc       write(*,*) 'lc',erx,ery,erz
cc       write(*,*) 'lc' ,rtt
       return
       end
       
       Subroutine Volcov(Rxx,Ryy,Rzz,Vxx,Vyy,Vzz,tt,dt,ww,mu,E0,kk)
       real*8 Rxx, Ryy, Rzz, Vxx, Vyy, Vzz
       real*8 E0(3), mu, ww, sins, coss
       real*8 tt, dt ,kk
cc       write(*,*) 'i am in Volcov',Rxx,Ryy,Rzz
cc       write(*,*) tt
       sins=Sin(ww*tt)-Sin(ww*(tt+dt))
       coss=Cos(ww*tt)-Cos(ww*(tt+dt))
cc       write(*,*) sins,coss,'after sins',E0(1),E0(2),E0(3)
       Rxx=Rxx-kk*(E0(1)/mu/ww**2.0)*coss+Vxx*dt
       Ryy=Ryy-kk*(E0(2)/mu/ww**2.0)*coss+Vyy*dt
       Rzz=Rzz-kk*(E0(3)/mu/ww**2.0)*coss+Vzz*dt
cc       write(*,*) Rxx,Ryy,Rzz,'before vel'
       Vxx=Vxx+kk*(E0(1)/mu/ww)*sins
       Vyy=Vyy+kk*(E0(2)/mu/ww)*sins
       Vzz=Vzz+kk*(E0(3)/mu/ww)*sins
cc       write(*,*) 'leaving Volcov'
       return
       end
       
       
       
       real function atanh(x)
       real*8 x
       atanh=0.5*log((x+1.0)/(1.0-x))
       return
       end
 

       subroutine  accels(xyze,xyzi,acce,acci,ionhist,ionchar,
     &  epairlist,
     &   ifeiint0,nepair,eiasm,nele,nion,outel,corechar)


!***calculates split operator corrected interation matrix******
!************************************************************** 
!************also finds a time step for next jump**************
       integer outel(512)
       integer nele
       integer nion
       integer eiasm(512)
       integer nepair
       integer epairlist(512,2)
       integer ionhist(512)
       integer    h
       integer     ninte
       integer     intel(512)
       real*8 acci(64,3)
       real*8 acce(512,3)
       real*8 veli(64,3)
       real*8 vele(512,3)
       real*8 xyzi(64,3)
       real*8 xyze(512,3)
       real*8 eeint(512,512)
       real*8 iiint(64,64)
       real*8 ieint(64,512)
       real*8 eiint(512,64)
       real*8 eeintfu(512,512)
       real*8 iiintfu(64,64)
       real*8 ieintfu(64,512)
       real*8 eiintfu(512,64)
       real*8 eedist(512,512)
       real*8 iidist(64,64)
       real*8 iedist(64,512)
       real*8 eidist(512,64)
       real*8  ifeeint(512,512)
       real*8 ifiiint(64,64)
       real*8 ifieint(64,512)
       real*8 ifeiint(512,64)
       real*8 ifeiint0(512,64)
       real*8 atomsize
       real*8 ionchar(64)
       real*8        zion
       real*8        dt
       real*8        mu
       real*8  corechar
       real*8 mi
       real*8 me
       me=1.0
       mi=20.1830*1836.0
       atomsize=2.5
cc       do i=1,nele
cc       write(*,*) ionchar(i)
cc       enddo
cc      stop

       j=0 
!       do i=1,512
       do i=1,nele
       if (ionhist(i).eq.1) then
       j=j+1
       intel(j)=i
       endif

       if ((ionhist(i).eq.0).and.(outel(i).eq.1)) then
       j=j+1
       intel(j)=i
       endif

       enddo
       ninte=j
cc       write(*,*) nele,'aaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaaa' 
!***********interaction matrix*********************************
!       do  h=1,512
!       do  n=1,512
!       do  m=1,512
!       do  l=1,512
!       do  k=1,512
!        do  n=1,1

!*******vectorization testing loop***************************** 
        
     
!*************direct Coulomb***********************************
!***********primary e-e interation matrix**********************
cc       write(*,*) ninte , 'ninte'
       do  i=1,ninte
       do  j=1,i
       eedist(intel(i),intel(j))=ABS((xyze(intel(i),1)-xyze(intel(j),1))
     & **2.0+
     &(xyze(intel(i),2)-xyze(intel(j),2))**2.0+(xyze(intel(i),3)-
     &xyze(intel(j),3))**2.0)**(1.0/2.0)
       eedist(intel(j),intel(i))=eedist(intel(i),intel(j))
       enddo
       enddo
       do  i=1,ninte
       do  j=1,i
       if (eedist(intel(i),intel(j)).gt.0.0) then
       eeint(intel(i),intel(j))=eedist(intel(i),intel(j))**(-3.0)
       endif 
       eeint(intel(j),intel(i))=eeint(intel(i),intel(j))
       enddo
       enddo
cc       write(*,*) ninte , 'ninte'
!***********primary i-i interaction matrix*********************
!      do  i=1,64
       do  i=1,nion
       do  j=1,i
       iidist(i,j)=ABS((xyzi(i,1)-xyzi(j,1))**2.0+(xyzi(i,2)-xyzi(j,2))
     &**2.0+(xyzi(i,3)-xyzi(j,3))**2.0)**(1.0/2.0)
       iidist(j,i)=iidist(i,j) 
       enddo
       enddo
       
!      do  i=1,64
       do  i=1,nion
       do  j=1,i
       if (iidist(i,j).gt.0.0) then
       iiint(i,j)=iidist(i,j)**(-3.0)
       endif
       iiint(j,i)=iiint(i,j)
       enddo
       enddo
c       write(*,*) iiint(1,1)
!***********primary i-e and e-i  interaction matrix*********************
!       do  i=1,64
       do  i=1,nion
!       do  j=1,512
       do  j=1,nele
       if((i.eq.eiasm(j)).or.(ionhist(j)+outel(j).gt.0)) then
!*****this is stupid for clearence and vec******************************
       iedist(i,j)=ABS((xyzi(i,1)-xyze(j,1))**
     &2.0+(xyzi(i,2)-xyze(j,2))**2.0+(xyzi(i,3)-xyze(j,3))
     &**2.0)**(1.0/2.0)
       endif
       enddo
       enddo
!       do  i=1,64
       do  i=1,nion
!       do  j=1,512
       do  j=1,nele
       eidist(j,i)=iedist(i,j)
       enddo
       enddo


!       do  i=1,64
       do  i=1,nion
!       do  j=1,512
       do  j=1,nele
       if((i.eq.eiasm(j)).or.(ionhist(j)+outel(j).gt.0)) then
       if (iedist(i,j).gt.0.0) then 
       ieint(i,j)=iedist(i,j)**(-3.0)
       endif
       if (eidist(j,i).gt.0.0) then 
       eiint(j,i)=eidist(j,i)**(-3.0)
       endif
       endif
       enddo
       enddo
c       write(*,*) ninte
!**************Condition determination**********************************
!***************Conditional interaction matrix**************************
!**********conditional e-e interaction matrix***************************
       do i=1,nele
       do j=1,nele
       ifeeint(i,j)=0.0
       enddo
       enddo
       do  i=1,ninte
       do  j=1,ninte
       ifeeint(intel(i),intel(j))=1.0
       enddo
       enddo
!**********conditional i-i interaction matrix***************************
!       do  i=1,64
       do i=1,nion
       do j=1,nion
       ifiiint(i,j)=0.0
       enddo
       enddo

       do  i=1,nion
       do  j=1,i
       if(i.ne.j) then
       ifiiint(i,j)=ionchar(i)*ionchar(j)
       ifiiint(j,i)=ifiiint(i,j)
       endif
       enddo
       enddo
!**********conditional i-e interaction matrix***************************
!       do  i=1,64
cc       write(3,*) 'next i-e'
       do  i=1,nion

       do  j=1,nele
       ifieint(i,j)=0.0
       enddo

       do  j=1,ninte
       ifieint(i,intel(j))=-ionchar(i)*1.0
       if(ionhist(intel(j)).eq.0) ifieint(i,intel(j))=0.0

       if((outel(intel(j)).eq.1).and.(i.ne.eiasm(intel(j))))
     & ifieint(i,intel(j))=-ionchar(i)*1.0
cc       write(3,*) ifieint(i,intel(j))
       enddo
       enddo
!*********conditional e-i interaction matrix****************************
!       do  i=1,512
       do  i=1,nele
!       do  j=1,64
       do  j=1,nion
       ifeiint(i,j)=-ifeiint0(i,j)*
     &(1-ionhist(i))-ionchar(j)*ionhist(i)
       if((outel(i).eq.1).and.(j.ne.eiasm(i)))
     & ifeiint(i,j)=-ionchar(j)
       enddo
       enddo
!***********************************************************************
!********calculates corrections to conditional interaction matrixes*****
!**********by  excluding interaction between paired particles***********
!************or makes swiss chees in interaction matrixes***************
!***********************************************************************
!*******************to all interaction matrixes*************************
!*********this is done to change easily to normal leap frog*************
       do i=1,nepair
       if((epairlist(i,1).lt.1000).and.(epairlist(i,2).lt.1000)) then
       ifeeint(epairlist(i,1),epairlist(i,2))=0.0
       ifeeint(epairlist(i,2),epairlist(i,1))=0.0
       endif
       
       if((epairlist(i,1).gt.1000).and.(epairlist(i,2).gt.1000)) then
       ifiiint(epairlist(i,1)-1000,epairlist(i,2)-1000)=0.0
       ifiiint(epairlist(i,2)-1000,epairlist(i,1)-1000)=0.0
       endif
       
       if((epairlist(i,1).gt.1000).and.(epairlist(i,2).lt.1000)) then
       ifieint(epairlist(i,1)-1000,epairlist(i,2))=0.0
       ifeiint(epairlist(i,2),epairlist(i,1)-1000)=0.0
       endif
       
       if((epairlist(i,1).lt.1000).and.(epairlist(i,2).gt.1000)) then
       ifeiint(epairlist(i,1),epairlist(i,2)-1000)=0.0
       ifieint(epairlist(i,2)-1000,epairlist(i,1))=0.0
       endif
       
       enddo

!*******************correction from assigned part***********************
!******************and from far electrons added*************************

       do i=1,nele
       do j=1,nion
       if ((eiasm(i).eq.j).and.(ionhist(i).eq.0)) ifeiint(i,j)=0.0 
       if (ionhist(i).eq.2) ifeiint(i,j)=0.0
       enddo
       enddo



!***********************************************************************
!******************Full interaction matrix******************************
!******************Full e-e interaction matrix**************************
       do  i=1,ninte
       do  j=1,ninte
       eeintfu(intel(i),intel(j))=eeint(intel(i),intel(j))*
     & ifeeint(intel(i),intel(j))
       enddo
       enddo
!******************Full i-i interaction matrix**************************
!       do  i=1,64
       do  i=1,nion
       do  j=1,i
       iiintfu(i,j)=iiint(i,j)*ifiiint(i,j) 
       iiintfu(j,i)=iiintfu(i,j)
       enddo
       enddo
!******************Full i-e interaction matrix**************************
!       do  i=1,64
       do  i=1,nion
       do  j=1,ninte
       ieintfu(i,intel(j))=ieint(i,intel(j))*ifieint(i,intel(j))
       enddo
       enddo
!******************Full e-i interaction matrix***************************
!       do  i=1,512
       do  i=1,nele
!       do  j=1,64
       do  j=1,nion
       eiintfu(i,j)=eiint(i,j)*ifeiint(i,j)
       enddo
       enddo
!******************Full e  acceleration matrix***************************
       do i=1,nele
       acce(i,1)=0.0
       acce(i,2)=0.0
       acce(i,3)=0.0
       enddo
cc       goto 10
!***************From e e Force*******************************************
       do i=1,ninte
       do k=1,ninte
       acce(intel(i),1)=acce(intel(i),1)+eeintfu(intel(i),intel(k))
     &*(xyze(intel(i),1)-xyze(intel(k),1))/me
       acce(intel(i),2)=acce(intel(i),2)+eeintfu(intel(i),intel(k))
     &*(xyze(intel(i),2)-xyze(intel(k),2))/me
       acce(intel(i),3)=acce(intel(i),3)+eeintfu(intel(i),intel(k))
     &*(xyze(intel(i),3)-xyze(intel(k),3))/me
       enddo
       enddo
!***************From e i  Force******************************************
!       do i=1,512
       do i=1,nele
!       do k=1,64
       do k=1,nion
       acce(i,1)=acce(i,1)+eiintfu(i,k)*(xyze(i,1)-xyzi(k,1))/me
       acce(i,2)=acce(i,2)+eiintfu(i,k)*(xyze(i,2)-xyzi(k,2))/me
       acce(i,3)=acce(i,3)+eiintfu(i,k)*(xyze(i,3)-xyzi(k,3))/me
       enddo
       enddo
cc10       continue
!******************Full i acceleration matrix****************************
       do i=1,nion
       acci(i,1)=0.0
       acci(i,2)=0.0
       acci(i,3)=0.0
       enddo
!***************From i e Force*******************************************
!       do i=1,64
       do i=1,nion
       do k=1,ninte
       acci(i,1)=acci(i,1)+ieintfu(i,intel(k))
     &*(xyzi(i,1)-xyze(intel(k),1))/mi
       acci(i,2)=acci(i,2)+ieintfu(i,intel(k))
     &*(xyzi(i,2)-xyze(intel(k),2))/mi
       acci(i,3)=acci(i,3)+ieintfu(i,intel(k))
     &*(xyzi(i,3)-xyze(intel(k),3))/mi
       enddo
       enddo
!***************From i i Force*******************************************
!       do i=1,64
       do i=1,nion
!       do k=1,64
       do k=1,nion
       acci(i,1)=acci(i,1)+iiintfu(i,k)*(xyzi(i,1)-xyzi(k,1))/mi
       acci(i,2)=acci(i,2)+iiintfu(i,k)*(xyzi(i,2)-xyzi(k,2))/mi
       acci(i,3)=acci(i,3)+iiintfu(i,k)*(xyzi(i,3)-xyzi(k,3))/mi
       enddo

cc       write(*,*) 'aciiiiiiiiiiiiiiiiiiiiiiiiiiiii'
cc       write(*,*) acci(i,1),acci(i,2),acci(i,3)
   

       enddo
cc10     continue 

!       enddo
!       enddo
!       write(*,*) h,h
!       enddo 
        return
       end 
       
       
       subroutine elli(ene,ee,xx,yy,zz,px,py,pz,kk,mm,sig)
       real*8 ene
       real*8 ee
       real*8 xx
       real*8 yy
       real*8 zz
       real*8 px 
       real*8 py
       real*8 pz
       real*8 mm
       real*8 sig
       real*8 r0
       real*8 ll
       real*8 kk
       r0=-kk/(2.0*ene)
       ll=(mm*kk*r0*(1.0-ee**2.0))**0.5*sig
       px=0.0
       py=ll/((1.0-ee)*r0)/mm
       pz=0.0
       xx=(1.0-ee)*r0
       yy=0.0
       zz=0.0
       write(*,*) xx,'xx'
       write(*,*) py,'py'
       write(*,*) py**2/2.0 -kk/xx
cc       stop
       return
       end 

       subroutine prepfull(eee,ener,kkk,nelat,nion,xyze,vele,veli,xyzi,
     & ionchar,ifeiint0,eiasm,outel,ionhist)
       real*8 xyzi(64,3)
       real*8 xyze(512,3) 
       real*8 vele(512,3)
       real*8 veli(64,3)
       real*8 eee(16)
       real*8 ener(16)
       real*8 kkk(16)
       real*8 axyz(16,3)
       real*8 avel(16,3) 
       real*8 xx
       real*8 yy
       real*8 zz
       real*8 px
       real*8 py
       real*8 pz
       real*8 kk
       real*8 me
       real*8 sig
       real*8 ee
       real*8 ene
       real*8 ionchar(64)
       real*8 ifeiint0(512,64)
       integer eiasm(512)
       integer nelat
       integer nion
       integer outel(512)
       integer ionhist(512)  
       me=1.0
       sig=+1.0
       do i=1,nelat
       ee=eee(i)
       ene=ener(i)
       kk=kkk(i)
       call elli(ene,ee,xx,yy,zz,px,py,pz,kk,me,sig)  
       axyz(i,1)=xx
       axyz(i,2)=yy
       axyz(i,3)=zz
       avel(i,1)=px
       avel(i,2)=py 
       avel(i,3)=pz 
       enddo 

       do i=1,nelat*nion
       do j=1, nion
       ifeiint0(i,j)=0.0
       enddo
       enddo

       do i=1,nion
       ionchar(i)=1.0
       do j=1,nelat
       k=nelat*(i-1)+j
       eiasm(k)=i
       write(*,*) eiasm(k),i,k
       ionhist(k)=0
       ifeiint0(k,i)=float(j)
       outel(k)=0
       if(j.eq.1) outel(k)=1
       xyze(k,1)=axyz(j,1)+xyzi(i,1)
       xyze(k,2)=axyz(j,2)+xyzi(i,2)
       xyze(k,3)=axyz(j,3)+xyzi(i,3)
       vele(k,1)=avel(j,1)+veli(i,1)
       vele(k,2)=avel(j,2)+veli(i,2)
       vele(k,3)=avel(j,3)+veli(i,3)
       enddo
       enddo
cc       stop
       return
       end 

       subroutine xray(xyze,vele,veli,xyzi,xr,xrt,ifcol,nele,nion,
     &    ionhist)
       real*8 xyzi(64,3)
       real*8 xyze(512,3)
       real*8 vele(512,3)
       real*8 veli(64,3)
       real*8 xr(10,64)
       real*8 xrt(10)
       real*8 bb(10)
       real*8 zz(10)
       real*8 ec(10)
       real*8 time 
       real*8 xatomsize
       real*8 ekin
       real*8 Pi 
       integer ifcol(512,64)
       integer ifcoll(512,64)
       integer nele
       integer nion
       integer ionhist(512)
       Pi=3.141592653
       xatomsize=2.91
       zz(1)=2.0
       zz(2)=8.0
       ec(1)=31.8713
       ec(2)=1.65
       bb(1)=0.35
       bb(2)=0.25
       do i=1,64
       xr(1,i)=0.0
       xr(2,i)=0.0
       enddo
       do i=1,nion
       do j=1,nele
       ifcoll(j,i)=0
!***prevents double counting of collisions*************************
       if (ionhist(j).eq.1) then 
       dist=ABS((xyzi(i,1)-xyze(j,1))**
     & 2.0+(xyzi(i,2)-xyze(j,2))**2.0+(xyzi(i,3)-xyze(j,3))
     & **2.0)**(1.0/2.0)
       if (dist.lt.xatomsize) ifcoll(j,i)=1 
 
       if ((dist.lt.xatomsize).and.(ifcol(j,i).eq.0)) then 
       ekin=ABS((veli(i,1)-vele(j,1))**
     & 2.0+(veli(i,2)-vele(j,2))**2.0+(veli(i,3)-vele(j,3))
     & **2.0)**(1.0/2.0)/2.0

       if(ekin.gt.ec(1)) then
       ekin=ekin/ec(1)
       xr(1,i)=xr(1,i)+Pi*bb(1)*zz(1)*log(ekin)/ekin/ec(1)**2.0
       endif
       if(ekin.gt.ec(2)) then
       ekin=ekin/ec(2)
       xr(2,i)=xr(2,i)+Pi*bb(2)*zz(2)*log(ekin)/ekin/ec(2)**2.0
       endif
       xr(1,i)=xr(1,i)/Pi/xatomsize**2.0
       xr(2,i)=xr(2,i)/Pi/xatomsize**2.0
cc       write(4,*) '**********************'
       endif
       endif
cc       write(4,*) ifcoll(j,i),'**',dist
       enddo
       enddo

       xrt(1)=0.0
       xrt(2)=0.0
       do i=1,64
       xrt(1)=xrt(1)+xr(1,i)
       xrt(2)=xrt(2)+xr(2,i)
       enddo
 
       do i=1,nion
       do j=1,nele
       ifcol(j,i)=ifcoll(j,i)
       enddo
       enddo
cc       write(4,*) xrt(1),xrt(2)
       return
       end












  
