        Program main
        implicit none
        integer(4),parameter :: n=2     ! dimension of nuclear coordinates
        real(8)              :: mass    ! mass
        real(8)              :: delta   ! time step
        real(8)              :: afric   ! friction parameter a=alpha*delta/2
        real(8)              :: rp(n),r0(n),rm(n) ! positions
        real(8)              :: force(n)
        real(8)              :: etot,epot,ekin
        real(8)              :: time
        integer(4)           :: niter     ! number of iterations
        integer(4)           :: iter0
        integer(4)           :: nfilt=10  ! fortran unit for trajectory file
        integer(4)           :: nfils=11  ! fortran unir for state file
        character(64)        :: trajectoryfile
        character(64)        :: statefile
        real(8)              :: svar1,svar2,svar3
        integer(4)           :: iter,i
!       ************************************************************************
!
!       ========================================================================
!       == variables set here are overwritten by the content of the state file==
!       ========================================================================
        mass=1.d0
        delta=1.d-1
        afric=0.d-1
        niter=1000
        r0(:)=(/1.d0,2.d0/)
        rm(:)=r0(:)
        trajectoryfile='tra.out'
        statefile='state.out'
        iter0=0
!
!       ========================================================================
!       == read state file                                                    ==
!       ========================================================================
        read(*,*)mass,delta,afric,niter
        read(*,*)r0,rm,iter0
        read(*,*)trajectoryfile,statefile
!
!       ========================================================================
!       == start molecular dynamics loop                                      ==
!       ========================================================================
!       == connect trajectoryfile to fortran unit nfilt
        OPEN(NFILT,FILE=TRAJECTORYFILE,ACTION='WRITE',ACCESS='APPEND')
        do iter=1,niter
          time=delta*real(iter0+iter,kind=8)
!
!         == calculate potential energy and forces =============================
          call energyandforce(n,r0,epot,force)
!
!         == propagate: determine position of next time step ===================
          svar1=2.d0/(1.d0+afric)
          svar2=-(1.d0-afric)/(1.d0+afric)
          svar3=delta**2/(mass*(1.d0+afric))
          do i=1,n
            rp(i)=r0(i)*svar1+rm(i)*svar2+force(i)*svar3
          enddo
!
!         == kinetic energy ====================================================
          ekin=0.d0
          do i=1,n
            ekin=ekin+0.5d0*mass*( (rp(i)-rm(i))/(2.d0*delta) )**2       
          enddo
!
!         == print information =================================================
          write(*,fmt='(4f10.5)')time,ekin,epot,ekin+epot
          write(nfilt,*)time,r0
!
!         == switch ============================================================
          do i=1,n
            rm(i)=r0(i)
            r0(i)=rp(i)
          enddo
        enddo
        close(nfilt) ! disconnect trajectory file
!
!       ========================================================================
!       == write state file                                                   ==
!       ========================================================================
        open(nfils,file=statefile)
        write(nfils,*)mass,delta,afric,niter
        write(nfils,*)r0,rm,iter0+niter
        write(nfils,*)"'"//trim(trajectoryfile)//"' '"//trim(statefile)//"'" 
        close(nfils)
        stop
        end
!
!       ........................................................................
        subroutine energyandforce(n,r,epot,force)
!       ************************************************************************
!       ** calculates potential energy and forces for positions r             **
!       ************************************************************************
        implicit none
        real(8)   ,parameter  :: c=1.d0         ! force constant
        real(8)   ,parameter  :: gamma=1.d0     ! electron-phonon coupling
        real(8)   ,parameter  :: t=0.1d0        ! hopping parameter
        real(8)   ,parameter  :: eb=gamma**2/(2*c) ! polaron binding energy
        integer(4),intent(in) :: n              ! number of atomic positions 
        real(8)   ,intent(in) :: r(n)           ! atomic coordinates
        real(8)   ,intent(out):: epot           ! potential energy
        real(8)   ,intent(out):: force(n)       ! forces
        real(8)               :: s,q            ! symmetry-adapted coordinates
        real(8)               :: deds,dedq      ! derivative of Epot wrt s,q
        real(8)               :: dsdr1,dsdr2,dqdr1,dqdr2
!       ************************************************************************
!       == transform r to symmetry adapted coordinates (s,q) ===================
        s=c/gamma*(r(1)+r(2))-1.d0
        q=c/gamma*(r(1)-r(2))
!       == evaluate potential energy ===========================================
        epot=0.5d0*Eb*(q**2-1.d0)-sqrt((Eb*q)**2+t**2) + 0.5d0*eb*s**2
!       == derivatives wrt (s,q) ===============================================
        deds=eb*s
        dedq=eb*q-0.5d0/sqrt( (Eb*q)**2+t**2 )*2.d0*eb**2*q
!       == forces F=-depot/dr ==================================================
        dsdr1=c/gamma
        dsdr2=c/gamma
        dqdr1=c/gamma
        dqdr2=-c/gamma
        force(1)=-dedq*dqdr1-deds*dsdr1
        force(2)=-dedq*dqdr2-deds*dsdr2
        return
        end 
