#include "defs.h"

#define MINSTEPS_PCTRL 10


!***********************************************************************
!
! temperature control
!
!***********************************************************************

SUBROUTINE VerletTempControl(deltat_fs,istep,transat)

  USE typeparam
  USE VTPcontrol

  IMPLICIT NONE

  ! *** output variables in module VTPcontrol ***

  ! xi  = "thermal" friction coefficients

  ! *** input variables in module VTPcontrol ***

  ! VTPtime        = time steps for T-P program
  ! VTempCtrl      = choice of temperature control
  ! NHtemp         = desired temperature
  ! NHQ            = the parameters of the Nose-Hoover thermostat

  ! *** variables passed in or out ***

  ! transat     = kinetic energy per atom and atom type
  !               (used for temperature contro
  ! istep       = current MD step
  ! deltat_fs   = time step in fs
  INTEGER,INTENT(IN)   :: istep
  REAL*8,INTENT(IN)    :: transat(itypelow:itypehigh),deltat_fs

  ! *** variables for parallel operation ***

  ! myatoms = number of atoms that my node is responsible for
#include "para_common.h"
#include "casc_common.h"

  ! *** local variables ***
  ! transavg   = average kinetic energy between the current
  !              and the preceding step
  ! transold   = kinetic energy of the preceding configuration per atom type
  ! dK         = ratio of the two kinetic energies defined above minus 1
  REAL*8       :: transavg,dK
  INTEGER      :: I,J
  REAL*8,SAVE,ALLOCATABLE :: transold(:)

  ! the equations used herein have been implemented according to the rapid
  ! communication of Hoolian, De Groot, and Hoover, PRA 41(1990), 4552 

  IF (istep.eq.1) THEN
     DO I=itypelow,itypehigh
        xi(I)=0.0d0
     END DO
  END IF
     
  IF (.NOT.ALLOCATED(transold)) THEN
     allocate(transold(itypelow:itypehigh))
     DO I=itypelow,itypehigh
        transold(I)=transat(I)
     END DO
  END IF

  ! *** temperature program ***

  ! find right set of variables belonging to current time step
  J=TPprogtimesteps
  DO WHILE(istep.ge.VTPtime(J+1).AND.J.gt.1)
     J=J-1
  END DO
  
  ! *** Nose-Hoover thermostat ***
  
  DO I=itypelow,itypehigh
     IF (VTempCtrl(I,J).eq.1) THEN
        transavg=half*(transat(I)+transold(I))
        transold(I)=transat(I)
        dK=two/(three*NHtemp(I,J))*transavg-one
        xi(i)=xi(i)+four*dK*deltat_fs/NHQ(I,J)
     ELSE
        xi(i)=0
     END IF
  END DO
  
END SUBROUTINE VerletTempControl




!***********************************************************************
!
! integrate the volume
!
!***********************************************************************

SUBROUTINE VerletVolumeControl(deltat_fs,istep,v,v0)

  USE typeparam
  USE VTPcontrol

  IMPLICIT NONE

  ! *** output variables in module VTPcontrol ***

  ! NHvscale = scaling factor of box size

  ! *** input variables from module VTPcontrol ***

  ! VTPtime        = time steps for T-P program
  ! VPressCtrl     = choice of temperature and pressure control
  ! eps            = friction coefficients associated with pressure control

  ! *** variables passed in or out ***

  ! istep       = current MD step
  ! deltat_fs   = time step in fs
  ! v,v0        = volumes (cf. dh, dhprev in mdx.f90)
  !               v=v(t+dt), v0=v(t)
  INTEGER,INTENT(IN)   :: istep
  REAL*8,INTENT(IN)    :: deltat_fs
  REAL*8,INTENT(INOUT) :: v,v0

  ! *** variables for parallel operation ***

  ! myatoms = number of atoms that my node is responsible for
#include "para_common.h"
#include "casc_common.h"

  ! *** local variables ***
  ! vold  = volumes, vold=v(t-dt)
  REAL*8,SAVE  :: vold
  INTEGER      :: I,J,K

  ! the equations used herein have been implemented according to the rapid
  ! communication of Hoolian, De Groot, and Hoover, PRA 41(1990), 4552 

  ! *** Nose-Hoover pressure control ***

  vold=v0
  v0=v
  ! integrate the volume
  IF (istep.eq.1) THEN
     v=v0+three*eps(1)*v0*deltat_fs
  ELSE
     v=vold+six*eps(1)*deltat_fs*v0
  END IF

  NHvscale(1)=(v/v0)**(one/three)
  NHvscale(2)=NHvscale(1) ; NHvscale(3)=NHvscale(1)
           
END SUBROUTINE VerletVolumeControl



!***********************************************************************
!
! pressure control
!
!***********************************************************************
SUBROUTINE VerletPressControl(deltat_fs,istep,trans,Pxx,Pyy,Pzz,v0)

  USE typeparam
  USE VTPcontrol

  IMPLICIT NONE

  ! *** output variables in module VTPcontrol ***

  ! eps = friction coefficient associated with pressure control
  ! NHvscale = scaling factor of box size

  ! *** input variables in module VTPcontrol ***

  ! VTPtime       = time steps for T-P program
  ! VPressCtrl    = choice of temperature and pressure control
  ! NHpress       = desired temperature and pressure
  ! NHpress[XYZ]  = desired pressure in xyz directions
  ! NHtau         = the parameters of the Nose-Hoover T and P control

  ! *** variables passed in or out ***

  ! trans       = kinetic energy per atom (used for pressure control)
  ! P[xyz]      = pressure components in xyz directions [eV/A^3]
  ! istep       = current MD step
  ! deltat_fs   = time step in fs
  ! v0          = current volume (cf. dh, dhprev in mdx.f90)
  INTEGER,INTENT(IN)   :: istep
  REAL*8,INTENT(IN)    :: trans,Pxx,Pyy,Pzz,deltat_fs,v0

  ! *** variables for parallel operation ***

  ! myatoms = number of atoms that my node is responsible for
#include "para_common.h"
#include "casc_common.h"

  ! *** local variables ***

  ! dP          = difference between desired and actual pressure
  ! pressure    = the actual pressure
  ! tausqinv    = the inverse of the squared pressure control parameter tau
  ! eps0,epsold = strain rates, eps0=eps(t), epsold=eps(t-dt)
  !               eps from module VTPcontrol: eps=eps(t+dt)
  REAL*8       :: tausqinv,dP,pressure,eps0
  REAL*8,SAVE  :: epsold(1:3)
  INTEGER      :: I,J,K

  ! the equations used herein have been implemented according to the rapid
  ! communication of Holian, De Groot, and Hoover, PRA 41(1990), 4552 

  IF (istep.ge.MINSTEPS_PCTRL) THEN
  
     ! *** pressure program ***

     ! find right set of variables belonging to current time step
     J=TPprogtimesteps
     DO WHILE(istep.ge.VTPtime(J+1).AND.J.gt.1)
        J=J-1
     END DO

     ! *** Nose-Hoover pressure control ***

     ! isotropic (hydrostatic) pressure
     IF (VPressCtrl(J).eq.1) THEN
        tausqinv=one/(NHtau(J)*NHtau(J))

        ! calculate the pressure difference in units of eV
        pressure=(Pxx+Pyy+Pzz)/three
        dP=pressure-NHpress(J)
        
        ! integrate the strain rate
        epsold(1)=eps0
        eps0=eps(1)
        IF (istep.eq.MINSTEPS_PCTRL) THEN
           eps(1)=eps0+three/two*v0*dP*tausqinv/trans*deltat_fs
        ELSE
           eps(1)=epsold(1)+three*v0*dP*tausqinv*deltat_fs/trans
        END IF

        eps(2)=eps(1) ; eps(3)=eps(1)
        epsold(2)=epsold(1) ; epsold(3)=epsold(1)

     ! direction dependent pressure
     ELSE IF (VPressCtrl(J).eq.1) THEN
        tausqinv=one/(NHtau(J)*NHtau(J))

        ! *** Pxx ***
        ! calculate the pressure difference in units of eV
        dP=Pxx-NHpressX(J)
        ! integrate the strain rate
        epsold(1)=eps0
        eps0=eps(1)
        IF (istep.eq.MINSTEPS_PCTRL) THEN
           eps(1)=eps0+three/two*v0*dP*tausqinv/trans*deltat_fs
        ELSE
           eps(1)=epsold(1)+three*v0*dP*tausqinv*deltat_fs/trans
        END IF

        ! *** Pyy ***
        ! calculate the pressure difference in units of eV
        dP=Pyy-NHpressY(J)
        ! integrate the strain rate
        epsold(2)=eps0
        eps0=eps(2)
        IF (istep.eq.MINSTEPS_PCTRL) THEN
           eps(2)=eps0+three/two*v0*dP*tausqinv/trans*deltat_fs
        ELSE
           eps(2)=epsold(2)+three*v0*dP*tausqinv*deltat_fs/trans
        END IF

        ! *** Pxx ***
        ! calculate the pressure difference in units of eV
        dP=Pzz-NHpressZ(J)
        ! integrate the strain rate
        epsold(3)=eps0
        eps0=eps(3)
        IF (istep.eq.MINSTEPS_PCTRL) THEN
           eps(3)=eps0+three/two*v0*dP*tausqinv/trans*deltat_fs
        ELSE
           eps(3)=epsold(3)+three*v0*dP*tausqinv*deltat_fs/trans
        END IF

     ELSE
        eps(1)=0    ; eps(2)=0    ; eps(3)=0
        epsold(1)=0 ; epsold(2)=0 ; epsold(3)=0
     END IF
  END IF

END SUBROUTINE VerletPressControl


!***********************************************************************
!
! Integrate Equations of Motion with the Verlet Algorithm
! 
!***********************************************************************

SUBROUTINE VerletIntegr(x0,x1,x2,x0nei,xnp,atype,box,boxsq,pbc,    &
     xneimax,xneimax2,delta,deltas,firstrun,                       &
     mtemp,heat,heatbrdr,Fzmaxz,FzmaxYz,Ekin,Elosssum,             &
     trans,transat,transin,natomsin,transbrdr,nbrdr,transv,        &
     Emax,Emaxbrdr,Emaxnosput,vmax,vmaxnosput,                     &
     tscaleth,tscalzmin,tscalzmax,tscalxout,tscalyout,             &
     timesputlim,tscalsputlim,fdamp)
  
  USE typeparam
  USE VTPcontrol, ONLY: eps,xi
  
  IMPLICIT NONE
  
  ! *** variables passed in and out ***

  ! x0,x1,x2  = positions and its time derivatives
  ! x0nei     = positions of nearest neighbors
  ! Ekin      = kinetic energy per atom
  ! Elosssum  = kinetic energy per atom lost in quenching process
  REAL*8,INTENT(INOUT)  :: x0(*),x1(*),x2(*),x0nei(*)
  REAL*8,INTENT(INOUT)  :: Ekin(*),Elosssum(*)

  ! *** variables passed in ***

  ! atype        = type of atoms
  ! mtemp        = parameter for temperature control (Gear)
  ! temp         = parameter for temperature control (Verlet)
  ! xnp          = forces
  ! box          = size of box in Angstroem
  ! pbc          = periodic boundary conditions (0.0 or 1.0)
  ! Fzmaxz, FzmaxYz = related to external forces
  ! delta(i)     = delta_fs/timeunit(i)
  ! deltas(i)    = delta(i)*delta(i)
  ! heat         = parameter for quenching of velocities
  ! heatbrdr     = quenching of broder atoms
  ! tscaleth     = skin thickness for definition of "border"
  ! Emaxbrdr     = maximum energy for border atoms
  ! transbrdr    = kinetic energy of border atoms
  ! nbrdr        = number of border atoms
  ! tscalzmin,tscalzmax = to define "border"
  ! timesputlim
  ! tscalsputlim = thickness of a layer in sputter simulation (--> border)
  ! fdamp        = damping force (quenching of velocities)
  ! tscal[xy]out = another way to assign the "at-border attribute" to atoms
  ! fc           = atom type dependent friction coefficient
  !                (effect of temperature and pressure control)
  ! firstrun  = used to set x0old,x0neiold values to zero in the first run
  LOGICAL,INTENT(IN) :: firstrun
  INTEGER,INTENT(IN) :: atype(*),mtemp
  REAL*8,INTENT(IN)  :: xnp(*),pbc(3),box(3),boxsq(3),Fzmaxz,FzmaxYz
  REAL*8,INTENT(IN)  :: delta(itypelow:itypehigh),deltas(itypelow:itypehigh)
  REAL*8,INTENT(IN)  :: heat,heatbrdr,tscaleth,Emaxbrdr
  REAL*8,INTENT(IN)  :: tscalzmin,tscalzmax,timesputlim
  REAL*8,INTENT(IN)  :: tscalsputlim,fdamp      
  REAL*8,INTENT(IN)  :: tscalxout,tscalyout

  ! *** variables passed out ***

  ! trans    = kinetic energy
  ! transv   = kinetic energy per direction
  ! transin  = kinetic energy of atoms within a given box
  ! transat  = kinetic energy per atom type
  ! Emax,vmax... = maximum values
  ! nbrdr    = number of atoms at the border
  ! natomsin = number of atoms with a given box
  REAL*8,INTENT(OUT)   :: trans,transin,transv(3),transbrdr
  REAL*8,INTENT(OUT)   :: transat(itypelow:itypehigh)
  REAL*8,INTENT(OUT)   :: Emax,vmax,Emaxnosput,vmaxnosput
  REAL*8,INTENT(OUT)   :: xneimax,xneimax2
  INTEGER,INTENT(OUT)  :: nbrdr,natomsin

  ! *** variables for parallel operation ***

  ! myatoms = number of atoms that my node is responsible for
#include "para_common.h"
#include "casc_common.h"

  ! *** local variables and constants ***

  REAL*8       :: tlim(3),E,v,help1,help2,heatfact
  REAL*8       :: x0new(3),x1store(3),Ebefore,Eafter
  REAL*8,SAVE  :: x0old(3*NPMAX),x0neiold(3*NPMAX)
  REAL*8       :: deltasq,d,boxsqoverdeltasq(3),fcj(1:3)
  INTEGER      :: i,j,j2,typej,i1,i2,i3,nscale
  LOGICAL      :: ifixed,isborderat
  ! Variables for neighbor list displacement calc.
  REAL*8       :: xneisqmax,xneisqmax2,xneisq

  ! tracking the number of calls to this function
  INTEGER itime
  DATA itime /0/
  SAVE itime

  ! ------------------------------------------------------------------
  ! ------------------------------------------------------------------
  ! ------------------------------------------------------------------

  ! *** count number of times subroutine is called ***

  itime=itime+1  

  ! *** set initial values ***
  
  trans=zero
  transin=zero
  natomsin=0
  transv(1)=zero; transv(2)=zero; transv(3)=zero; 
  transbrdr=zero
  DO j=itypelow,itypehigh
     transat(j)=0.0d0
  END DO   

  nscale=0
  nbrdr=0
  DO i=1,3
     ! thickness of border layer
     tlim(i)=half-tscaleth/box(i)
  ENDDO
  
  ! maximum kinetic energies
  Emax=-1.0d0
  Emaxnosput=-1.0d0
  ! maximum velocities
  vmax=-1.0d0
  vmaxnosput=-1.0d0
  ! max and next to max displacement
  xneisqmax=0.0d0
  xneisqmax2=0.0d0
  ! quenching velocities
  heatfact=1.0-heat*heat

  ! set initial values for 'old' arrays
  IF (firstrun) THEN
     DO i=1,3*NPMAX
       x0old(i)=x0(i)
       x0neiold(i)=x0nei(i)
     ENDDO
  ENDIF

  AL_DEBUG('x0 sum before Verlet integration',x0,3*myatoms)
  AL_DEBUG('xnp sum before Verlet integration',xnp,3*myatoms)
  AL_DEBUG('x2 sum before Verlet integration',x2,3*myatoms)

  ! *** loop over atoms on this processor ***

  DO 250 j=1,myatoms

     ! *** quantities depending on the type of the current atom ***

     v=zero
     typej=ABS(atype(j))
     ! time steps
     d=delta(typej)
     deltasq=d*d
     fcj(1)=d*(xi(typej)+eps(1))
     fcj(2)=d*(xi(typej)+eps(2))
     fcj(3)=d*(xi(typej)+eps(3))
     xneisq=zero
     ! scaled boxsizes
     boxsqoverdeltasq(1)=boxsq(1)/deltasq
     boxsqoverdeltasq(2)=boxsq(2)/deltasq
     boxsqoverdeltasq(3)=boxsq(3)/deltasq        
     ! indices for coordinates
     i1=j*3-2
     i2=j*3-1
     i3=j*3-0
     ! free or fixed atom
     ifixed=.false.
     if (atype(j)<0) ifixed=.true.

     ! *** KN: Determine what are border atoms ***
          
     ! This is terribly complex - but there is a reason for everything
     
     ! initial "assumption"
     isborderat=.false.     
     
     ! part 1:  arguments for being a border atom
     if (x0(i3)*box(3)>tscalzmin.and.x0(i3)*box(3)<tscalzmax) then
        isborderat=.true.
     else if (x0(i1)*box(1)<-tscalxout .or. x0(i1)*box(1)>tscalxout) then
        isborderat=.true.
     else if (x0(i2)*box(2)<-tscalyout .or. x0(i2)*box(2)>tscalyout) then
        isborderat=.true.
     else if (x0(i1)>tlim(1).or.x0(i1)<=-tlim(1) .or. &
          x0(i2)>tlim(2).or.x0(i2)<=-tlim(2) .or.     &
          x0(i3)>tlim(3).or.x0(i3)<=-tlim(3)) then
        ! don't scale at open borders
        if ((mtemp==7.or.mtemp==8).and. pbc(1)/=1.0 .and.         &
             (x0(i1)>tlim(1).or.x0(i1)<-tlim(1))) goto 20
        if ((mtemp==7.or.mtemp==8).and.pbc(2)/= 1.0 .and.         &
             (x0(i2)>tlim(2).or.x0(i2)<-tlim(2))) goto 20
        if ((mtemp==7.or.mtemp==8).and.pbc(3)/= 1.0 .and.         &
             (x0(i3)>tlim(3).or.x0(i3)<-tlim(3))) goto 20
        
        isborderat=.true.
     endif
20   continue
     
     ! part 2: arguments for not being "regarded" as a border atom
     if (ifixed) isborderat=.false.
     if (myproc == irecproc .and. j == irec) isborderat=.false.
     if (x0(i1) < -tscalsputlim .or. x0(i1) > tscalsputlim .or.      &
          x0(i2) < -tscalsputlim .or. x0(i2) > tscalsputlim .or.     &
          x0(i3) < -tscalsputlim .or. x0(i3) > tscalsputlim ) then
        isborderat=.false.
     endif
     
     ! *** temperature control for border atoms ***

     ! Now, that we know the status of our comrade, we can
     ! start playing the thermostat game but before we would
     ! like to know the energy so that we can find the difference
     ! due to the thermostat later
     help2=      x1(i1)*x1(i1)*boxsqoverdeltasq(1)
     help2=help2+x1(i2)*x1(i2)*boxsqoverdeltasq(2)
     help2=help2+x1(i3)*x1(i3)*boxsqoverdeltasq(3)
     Ebefore=half*help2
     
     if (isborderat) then
        ! heatbrdr is = 1 except if mtemp=5, 7 or 8,
        ! as determined in Temp_Control().
        x1(i1)=heatbrdr*x1(i1)
        x1(i2)=heatbrdr*x1(i2)
        x1(i3)=heatbrdr*x1(i3)
        nscale=nscale+1
     endif
     
     ! now, calculate energy after border temperature scale
     help2=      x1(i1)*x1(i1)*boxsqoverdeltasq(1)
     help2=help2+x1(i2)*x1(i2)*boxsqoverdeltasq(2)
     help2=help2+x1(i3)*x1(i3)*boxsqoverdeltasq(3)
     Eafter=half*help2;
     
     ! get heat loss sums
     if (isborderat) then
        Elosssum(1)=Elosssum(1)+Ebefore-Eafter
     endif
     
     ! now, use Eafter as start value of energy calculation
     ! for overall energy loss (i.e. including non-border atoms)
     Ebefore=Eafter
        
     ! *** loop over directions ***
     
     DO 200 j2=1,3
        ! position in x0,x1... arrays
        i=(j-1)*3+j2
        
        ! quenching
        x1(i)=heat*x1(i)
        ! store this velocity for temperature-scaling-energy-loss calculation
        x1store(j2)=x1(i)

        ! *** integrate equations of motion ***

        ! Only place where time step enters solution of equations of motion!!
        x2(i)=two*deltas(typej)*xnp(i)
        ! the friction coefficient as calculated in the temperature
        ! and pressure control
        x2(i)=x2(i)-fcj(j2)*x1(i)
     
        ! Verlet algorithm
        IF (.NOT. ifixed .AND..NOT.(FzmaxYz/=zero &
             .AND. j2==3 .AND. x0(i)*box(3)>=Fzmaxz)) THEN
           if (firstrun) THEN 
              ! integrate positions
              x0new(j2)=x0(i)+x1(i)+x2(i)
           ELSE

              IF (x2(i).eq.zero) THEN
                 ! integrate positions
                 x0new(j2)=x0(i)+x1(i)
              ELSE
                 ! integrate positions
                 x0new(j2)=two*x0(i)-x0old(i)+x2(i)
                 ! integrate velocities: central difference
                 x1(i)=half*(x0new(j2)-x0old(i))      
              END IF
           END IF
           x0old(i)=x0(i)
           x0(i)=x0new(j2)
           xneisq=xneisq+x0nei(i)*x0nei(i)*boxsq(j2)
           ! integrate positions in neighbor list, too
           x0new(j2)=two*x0nei(i)-x0neiold(i)+x2(i)
           x0neiold(i)=x0nei(i)
           x0nei(i)=x0new(j2)
        ELSE
           x1(i)=zero
        ENDIF

        ! *** thermostats ***

        ! thermostat type 2:
        ! quench - if energy goes up, set all velocities to zero
        IF (mtemp==2) x1(i)=zero

        ! thermostat type 9:
        ! not documented in README.DOCUMENTATION
        if (mtemp==9.and.j2==1) then
           if (x1(i)*x2(i)+x1(i+1)*x2(i+1)+x1(i+2)*x2(i+2)<zero) then
              x1(i  )=fdamp*x1(i  )
              x1(i+1)=fdamp*x1(i+1)
              x1(i+2)=fdamp*x1(i+2)
           end if
        end if

        ! part of calculation of kinetic energies
        help2=x1(i)*x1(i)*boxsqoverdeltasq(j2) ! boxs(j2)/deltasq

        ! *** periodic boundaries ***
        
        IF (pbc(j2) .EQ. 1.0 .AND. ABS(x0(i)) .GE. 1.5d0) THEN
           WRITE(6,100) myproc,i,x0(i),xnp(i),x2(i)
           WRITE(7,100) myproc,i,x0(i),xnp(i),x2(i)
100        FORMAT('PROBLEM IN VERLET: node,i,x0(i),xnp(i),x2(i)=',  &
                2i5,3f20.10)
        ENDIF
        
        IF (x0(i) .GE. half) THEN
           x0(i)=x0(i)-pbc(j2)
        ELSEIF (x0(i) < -half) THEN
           x0(i)=x0(i)+pbc(j2)
        ENDIF

        ! adding up kinetic contributions to the energy per direction
        ! used for the calculation of internal pressure
        transv(j2)=transv(j2)+help2  
        v=v+help2                    ! calc of total

200  END DO  ! End of loop over dimensions 1 - 3

     ! calculate energy after temperature scale, but excluding the 
     ! solution of the equations of motion
     help2=      x1store(1)*x1store(1)*boxsqoverdeltasq(1)
     help2=help2+x1store(2)*x1store(2)*boxsqoverdeltasq(2)
     help2=help2+x1store(3)*x1store(3)*boxsqoverdeltasq(3)
     Eafter=half*help2;

     ! Get heat loss sum 
     Elosssum(2)=Elosssum(2)+Ebefore-Eafter
	 
     ! *** KE: get max and next to max displacement ***
     
     if (xneisq > xneisqmax) then
        xneisqmax2=xneisqmax
        xneisqmax=xneisq
        ! print *,'max1',xneisq,xneisqmax,xneisqmax2  
     else 
        if (xneisq > xneisqmax2) then
           xneisqmax2=xneisq
           ! print *,'max2',xneisq,xneisqmax,xneisqmax2
        endif
     endif
	 
     ! *** kinetic energies ***

     E=half*v
     v=sqrt(v)*vunit(typej)
     Ekin(j)=E
     trans=trans+E
     transat(typej)=transat(typej)+E

     ! *** KE: get max E and v (be careful with the units !) ***
	 
     if (E > Emax) Emax=E
     if (v > vmax) vmax=v
     
     if (.not.(x0(i1)<-timesputlim .or. x0(i1)>timesputlim .or.        &
          x0(i2)<-timesputlim .or. x0(i2)>timesputlim .or.             &
          x0(i3)<-timesputlim .or. x0(i3)>timesputlim)) then
        transin=transin+E
        natomsin=natomsin+1
     endif
	    
     if (.not.(x0(i1)<-timesputlim .or. x0(i1)>timesputlim .or.    &
          x0(i2)<-timesputlim .or. x0(i2)>timesputlim .or.         &
          x0(i3)<-timesputlim .or. x0(i3)>timesputlim)) then
        if (E > Emaxnosput) Emaxnosput=E
        if (v > vmaxnosput) vmaxnosput=v
     endif


     ! *** KN: Determine what are border atoms ***

     ! These criteria should match the ones above exactly.

     ! initial "assumption"
     isborderat=.false.     

     ! part 1:  arguments for being a border atom
     if (x0(i3)*box(3)>tscalzmin.and.x0(i3)*box(3)<tscalzmax) then
        isborderat=.true.
     else if (x0(i1)*box(1)<-tscalxout .or. x0(i1)*box(1)>tscalxout) then
        isborderat=.true.
     else if (x0(i2)*box(2)<-tscalyout .or. x0(i2)*box(2)>tscalyout) then
        isborderat=.true.
     else if (x0(i1)>tlim(1).or.x0(i1)<=-tlim(1) .or.           &
          x0(i2)>tlim(2).or.x0(i2)<=-tlim(2) .or.               &
          x0(i3)>tlim(3).or.x0(i3)<=-tlim(3)) then
        ! Don't scale at open borders
        if ((mtemp==7.or.mtemp==8).and. pbc(1)/=1.0 .and.       &
             (x0(i1)>tlim(1).or.x0(i1)<-tlim(1))) goto 30
        if ((mtemp==7.or.mtemp==8).and.pbc(2)/= 1.0 .and.       &
             (x0(i2)>tlim(2).or.x0(i2)<-tlim(2))) goto 30
        if ((mtemp==7.or.mtemp==8).and.pbc(3)/= 1.0 .and.       &
             (x0(i3)>tlim(3).or.x0(i3)<-tlim(3))) goto 30
        
        isborderat=.true.
     endif
30   continue
     
     ! part 2: arguments for not being 'regarded' as a border atom
     if (ifixed) isborderat=.false.
     if (myproc == irecproc .and. j == irec) isborderat=.false.
     if (x0(i1) < -tscalsputlim .or. x0(i1) > tscalsputlim .or.      &
          x0(i2) < -tscalsputlim .or. x0(i2) > tscalsputlim .or.     &
          x0(i3) < -tscalsputlim .or. x0(i3) > tscalsputlim ) then
        isborderat=.false.
     endif
     
     ! *** get border kinetic energy ***
     
     if (isborderat) then
        transbrdr=transbrdr+E
        nbrdr=nbrdr+1
        if (E > Emaxbrdr) then
           WRITE(6,'(//,A,I8,A,F10.2,A,F10.2)') 'ERROR: atom # ',j,' Ekin ',E,' > Emaxbrdr',Emaxbrdr
           WRITE(6,'(A,3F12.3)') 'Atom position ',  &
                x0(i1)*box(1),x0(i2)*box(2),x0(i3)*box(3)
           MYNFLUSH(6)	       
           AL_STOP('Emaxbrdr exceeded',myproc)
        endif
     endif     

250 END DO     

  AL_DMAX(vmax,1)
  AL_DMAX(Emax,1)
  AL_DMAX(Emaxnosput,1)
  AL_DMAX(vmaxnosput,1)
  DPRINT *,'Emax',Emax,Emaxnosput
  
  xneimax=sqrt(xneisqmax)
  xneimax2=sqrt(xneisqmax2)
  
  AL_DMAX(xneimax,1)
  AL_DMAX(xneimax2,1)
  
  AL_DEBUG('x0 sum in Verlet is',x0,3*myatoms)

END SUBROUTINE VerletIntegr



!***********************************************************************
!
! calculate kinetic energy
!
!***********************************************************************

SUBROUTINE CalcKinEnergies(x1,atype,boxsq,delta,Ekin,trans,transv,transat)
  
  USE typeparam

  IMPLICIT NONE

  ! *** variables passed in and out ***
  
  ! x1...  = velocities
  ! boxsq  = squared size of box in Angstroem
  ! Ekin   = kinetic energies of atoms
  ! trans  = summed up kinetic energy
  ! transv = summed up kinetic energy contributions per coordinate
  ! transat = summed up kinetic energy per atom type
  REAL*8,INTENT(IN)  :: x1(*),boxsq(3),delta(itypelow:itypehigh)
  INTEGER,INTENT(IN) :: atype(*)
  REAL*8,INTENT(OUT) :: Ekin(*),trans,transv(3),transat(itypelow:itypehigh)

  ! *** variables for parallel operation ***
  
  ! myatoms = number of atoms that my node is responsible for
#include "para_common.h"
#include "casc_common.h"

  ! *** local variables and constants ***
  
  REAL*8  :: help1,help2,deltasq,d,boxsqoverdeltasq(3)
  INTEGER :: i,j,j2,typej
  
  ! ------------------------------------------------------------------
  ! ------------------------------------------------------------------
  ! ------------------------------------------------------------------

  ! *** loop over all atoms on this processor ***
  
  trans=zero
  transv(1)=0 ; transv(2)=0 ; transv(3)=0
  DO j=itypelow,itypehigh
     transat(j)=0.0d0
  END DO   

  DO j=1,myatoms
     typej=atype(j)
     d=delta(typej)
     deltasq=d*d
     help2=0
     DO j2=1,3
        boxsqoverdeltasq(j2)=boxsq(j2)/deltasq
        i=(j-1)*3+j2
        ! for calc of kinetic energy per atom
        help1=x1(i)*x1(i)*boxsqoverdeltasq(j2)
        help2=help2+help1
        ! add up kinetic contributions to the total energy per direction
        ! used for the calculation of internal pressure
        transv(j2)=transv(j2)+help1
     END DO
     Ekin(j)=half*help2
     trans=trans+Ekin(j)
     transat(typej)=transat(typej)+Ekin(j)
  END DO
     
END SUBROUTINE CalcKinEnergies
