!********************************************************************
! file:        meso.f90
! language:    f90
! author:      Emma Falck
! updated:     April, 23
! description: A program for the implementation of the
!              mesoscopic simulation algorithm described in
!              Malevanets et al., J. Chem. Phys. 112, 7260 (2000).
! version:     One polymer (NSE monomers) and NST solvent 
!              particles in 2D.
!********************************************************************
! The author of this software is E. Falck.  Copyright (c) 2002.
! Permission to use, copy, modify, and distribute this software for 
! any purpose without fee is hereby granted, provided that this entire
! notice is included in all copies of any software which is or includes
! a copy or modification of this software and in all copies of the 
! supporting documentation for such software.
! THIS SOFTWARE IS BEING PROVIDED "AS IS", WITHOUT ANY EXPRESS OR 
! IMPLIED WARRANTY. IN PARTICULAR, NEITHER THE AUTHOR NOR SOFTSIMU 
! MAKE ANY REPRESENTATION OR WARRANTY OF ANY KIND CONCERNING THE 
! MERCHANTABILITY OF THIS SOFTWARE OR ITS FITNESS FOR ANY PARTICULAR 
! PURPOSE.
!********************************************************************

PROGRAM mesosim

!--------------------------------------------------------------------
! The constants and parameters used in this program are 
! imported from the module "parameters". 
!--------------------------------------------------------------------

  USE parameters

  IMPLICIT NONE

!--------------------------------------------------------------------
! The following functions/subroutines are imported for random 
! number generation.
!--------------------------------------------------------------------

  EXTERNAL sgrnd
  REAL(DP), EXTERNAL :: grnd

!--------------------------------------------------------------------
! Variable declarations.
!--------------------------------------------------------------------

  INTEGER :: &
      checkstep, &   ! step index for checking quantities such as \beta
      colstep, &     ! collision step index
      mdstep, &      ! MD step index
      straccount, &  ! counter for stress-stress a.c. function
      tstrac, &      ! time index for stress-stress a.c. function
      tvcmac, &      ! time index for c. of m. velocity a.c. function
      vcmaccount     ! counter for c. of m. velocity a.c. function

  INTEGER, DIMENSION(1:NSE + 1) :: &
      point         ! index to the Verlet neighbor list
    
  INTEGER, DIMENSION(1:MAXNAB) :: &
      list          ! Verlet neighbour list

  REAL(DP) :: &
      pot, &         ! total potential energy of the system
      vir            ! virial for the system

  REAL(DP), DIMENSION(1:DIM) :: &
      centmass, &    ! center of mass of the polymer
      cmveloc        ! center of mass velocity 

  REAL(DP), DIMENSION(1:MAXCHECK) :: &
      poten, &       ! total potential energy of the system
      testbeta, &    ! computed values of 1 / (k_B T)
      virial         ! virial for the system

  REAL(DP), DIMENSION(1:COLSTEPS) :: &
      endtoendsq, &  ! end-to-end distance of the polymer
      rgs            ! radius of gyration squared

  INTEGER, DIMENSION(1:MAXSTRAC + 1) :: &
      tstracbound    ! for calculating stress-stress a.c. f. "on the fly"
  
  INTEGER, DIMENSION(-MAXSTRAC:MAXSTRAC) :: &
      stracloop      ! for calculating stress-stress a.c. f. "on the fly"

  REAL(DP), DIMENSION(0:MAXSTRAC - 1) :: &
      strac          ! stress-stress autocorrelation function (Malevanets) 

  REAL(DP), DIMENSION(0:MAXSTRAC - 2) :: &
      visc           ! viscosity (Malevanets)

  REAL(DP), DIMENSION(1:MAXSTRAC) :: &
      stracstore     ! stores MAXSTRAC stress-stress elements 

  INTEGER, DIMENSION(1:MAXVCMAC + 1) :: &
      tvcmacbound    ! for calculating c. of m. velocity a.c. f. "on the fly"
  
  INTEGER, DIMENSION(-MAXVCMAC:MAXVCMAC) :: &
      vcmacloop      ! for calculating c. of m. velocity a.c. f. "on the fly"

  REAL(DP), DIMENSION(0:MAXVCMAC - 1) :: &
      vcmac          ! center of mass velocity autocorrelation function 

  REAL(DP), DIMENSION(1:MAXVCMAC,1:DIM) :: &
      vcmacstore     ! stores MAXVCMAC center of mass velocities

  REAL(DP), DIMENSION(1:MAXCHECK,1:DIM) :: &
      vcm            ! center of mass velocity

  CHARACTER(LEN = MAXSTR) :: &
      filename, &    ! used in naming files
      tempname       ! used in naming files

!--------------------------------------------------------------------
! Quantities related to the collision volumes.
!--------------------------------------------------------------------

  REAL(DP), DIMENSION(1:NWS) :: &
      cvamount       ! amount of particles in collision volume

  REAL(DP), DIMENSION(1:NWS,1:DIM) :: &
      cvcoord, &     ! coordinates of collision volumes
      cvaveveloc     ! average velocities of collision volumes

  REAL(DP), DIMENSION(1:NWS,1:DIM,1:DIM) :: &
      cvrotation     ! rotation matrices for collision volumes

!--------------------------------------------------------------------
! Quantities related to solvent and solute particles.
!--------------------------------------------------------------------

  INTEGER, DIMENSION(1:NST) :: &
      cvind          ! indicates to which coll. vol. a particle belongs
  
  REAL(DP), DIMENSION(1:NST) :: &
      deltaxii       ! \xi_iy(\tau + 1) - \xi_iy(\tau)

  REAL(DP), DIMENSION(1:N,1:DIM) :: &
      acc, &         ! particle accelerations
      oldposit, &    ! particle positions
      posit, &       ! particle positions
      veloc          ! particle velocities

  REAL(DP), DIMENSION(1:NSE,1:DIM) :: &
      positnpbc      ! solute posit. without periodic boundary conditions

!--------------------------------------------------------------------
! The random number generator is initialized.
!--------------------------------------------------------------------

  CALL sgrnd(RSEED)

!--------------------------------------------------------------------
! The 2D collision volume coordinates are initialized.
!--------------------------------------------------------------------

  CALL init_cell_coords(cvcoord)

!--------------------------------------------------------------------
! The positions of the solvent and solute particles are initialized.
!--------------------------------------------------------------------

  CALL init_posit_se(posit,positnpbc)
  CALL init_posit_st(posit)

!--------------------------------------------------------------------
! The velocities of all particles are initialized.
!--------------------------------------------------------------------

  CALL init_veloc(veloc)

!--------------------------------------------------------------------
! The quantities needed in the computation of the stress-stress 
! autocorrelation function are initialized.
!--------------------------------------------------------------------

  CALL init_ac(1,MAXSTRAC,tstrac,straccount,tstracbound, &
      stracloop,strac)

!-------------------------------------------------------------------
! The quantities needed in the computation of the
! center of mass velocity autocorrelation function 
! of the polymer are initialized.
!-------------------------------------------------------------------

  CALL init_ac(1,MAXVCMAC,tvcmac,vcmaccount,tvcmacbound, &
      vcmacloop,vcmac)

!--------------------------------------------------------------------
! Additional (simple) initializations are done.
!--------------------------------------------------------------------

  checkstep = 1

!--------------------------------------------------------------------
! The system is brought to the desired temperature. 
!--------------------------------------------------------------------  

  DO colstep = 1, TEMPSTEPS

!--------------------------------------------------------------------
! The MD steps are taken. First we need to compute the "initial"
! accelerations. 
!--------------------------------------------------------------------

    CALL comp_acc(.true.,point,list,pot,vir,acc,posit)

    DO mdstep = 1, MDSTEPS

      CALL do_md(point,list,pot,vir,acc,posit,veloc,positnpbc)
      CALL scale_veloc(veloc)

    END DO

!--------------------------------------------------------------------
! The particles are folded back to the original simulation box.
!--------------------------------------------------------------------

    CALL fold_back(posit)

!--------------------------------------------------------------------
! The collision volume indices are updated. 
!--------------------------------------------------------------------

    CALL assign_cell_indices(cvind,posit)

!--------------------------------------------------------------------
! The collision step for the solvent is taken.
!--------------------------------------------------------------------

    CALL collide_st(cvamount,cvaveveloc,cvrotation,cvind,veloc)

  END DO

!--------------------------------------------------------------------
! The system is equilibrated. 
!--------------------------------------------------------------------  

  DO colstep = 1, EQSTEPS

!--------------------------------------------------------------------
! The MD steps are taken. 
!--------------------------------------------------------------------

    CALL comp_acc(.true.,point,list,pot,vir,acc,posit)

    DO mdstep = 1, MDSTEPS

      CALL do_md(point,list,pot,vir,acc,posit,veloc,positnpbc)

    END DO

!--------------------------------------------------------------------
! The particles are folded back to the original simulation box.
!--------------------------------------------------------------------

    CALL fold_back(posit)

!--------------------------------------------------------------------
! The collision volume indices are updated. 
!--------------------------------------------------------------------

    CALL assign_cell_indices(cvind,posit)

!--------------------------------------------------------------------
! The collision step for the solvent is taken.
!--------------------------------------------------------------------

    CALL collide_st(cvamount,cvaveveloc,cvrotation,cvind,veloc)

  END DO

!--------------------------------------------------------------------
! The actual simulation starts here. The dynamics is the same
! as above. In addition, relevant quantities are computed 
! and recorded.
!--------------------------------------------------------------------
 
  DO colstep = 1, COLSTEPS

!--------------------------------------------------------------------
! The configuration is saved every EQCONFSTEPS.
!--------------------------------------------------------------------

    IF (mod(colstep,CONFSTEPS) == 0) THEN

      write(tempname,NAMEFORM) colstep
      filename = trim(MPREFIX) // trim(adjustl(tempname(1:MAXSTR))) &
        // trim(POSTFIX)
      CALL output_configuration(filename,CONFFORM,positnpbc)

    END IF

    CALL comp_acc(.true.,point,list,pot,vir,acc,posit)

    DO mdstep = 1, MDSTEPS

      CALL do_md(point,list,pot,vir,acc,posit,veloc,positnpbc)

    END DO

!--------------------------------------------------------------------
! The displacement in the y direction from the previous collision 
! volume to the current one is computed. Notice that this must be
! done before the foldback.
!--------------------------------------------------------------------

    CALL comp_deltaxii(cvcoord,cvind,deltaxii,posit)

    CALL fold_back(posit)

    CALL assign_cell_indices(cvind,posit)

!--------------------------------------------------------------------
! The end-to-end distance is computed.
!--------------------------------------------------------------------

    CALL comp_endtoend(endtoendsq(colstep),positnpbc)

!--------------------------------------------------------------------
! The center of mass and the radius of gyration (squared) 
! of the polymer are computed. 
!--------------------------------------------------------------------

    CALL comp_centmass(centmass,positnpbc)
    CALL comp_radgyr(rgs(colstep),centmass,positnpbc)
      
!--------------------------------------------------------------------
! The quantities needed for the calculation of the 
! stress-stress autocorrelation function are computed and recorded.
!--------------------------------------------------------------------

    CALL comp_stress(stracstore(tstrac),deltaxii,veloc)

    CALL comp_ac(1,MAXSTRAC,1,tstrac,straccount, &
        tstracbound,stracloop,strac,stracstore)

!--------------------------------------------------------------------
! The quantities needed for the computation of the center of mass 
! velocity autocorrelation function of the polymer are taken 
! care of.
!--------------------------------------------------------------------

    CALL comp_cmveloc(cmveloc,veloc)
    vcmacstore(tvcmac,:) = cmveloc
    CALL comp_ac(1,MAXVCMAC,DIM,tvcmac,vcmaccount, &
        tvcmacbound,vcmacloop,vcmac,vcmacstore)

    CALL collide_st(cvamount,cvaveveloc,cvrotation,cvind,veloc)

!--------------------------------------------------------------------
! At certain intervals we store the potential, the virial, the
! temperature, and the center of mass velocity.
!--------------------------------------------------------------------

    IF (mod(colstep,CHECKSTEPS) == 0) THEN

      poten(checkstep) = pot
      virial(checkstep) = vir
      vcm(checkstep,:) = cmveloc

      CALL check_temp(testbeta(checkstep),veloc)

      checkstep = checkstep + 1

    END IF

  END DO

!--------------------------------------------------------------------
! The computation of the stress-stress autocorrelation function is
! completed. According to Malevanets et al., the viscosity is 
! computed from the stress-stress autocorrelation function.
!--------------------------------------------------------------------
  
  CALL fin_ac(1,MAXSTRAC,straccount,strac)
  CALL comp_visc(strac,visc)

!--------------------------------------------------------------------
! The computation of the center of mass velocity autocorrelation 
! function of the polymer is completed. Note that velocity
! autocorrelation functions are often scaled by the dimension.
!--------------------------------------------------------------------
  
  CALL fin_ac(1,MAXVCMAC,vcmaccount,vcmac)
  vcmac = vcmac / DIM 

!--------------------------------------------------------------------
! The output is written.
!--------------------------------------------------------------------

  CALL output_general('general.txt')
  CALL output_real_data('poten.dat',MAXCHECK,CHECKSTEPS * TAU, &
      CHECKSTEPS * TAU,poten)
  CALL output_real_data('virial.dat',MAXCHECK,CHECKSTEPS * TAU, &
      CHECKSTEPS * TAU,virial)
  CALL output_real_data('endtoendsq.dat',COLSTEPS,TAU,TAU,endtoendsq)
  CALL output_real_data('rgs.dat',COLSTEPS,TAU,TAU,rgs)
  CALL output_real_data('stressac.dat',MAXSTRAC,0._dp,TAU,strac)
  CALL output_real_data('viscosity.dat',MAXSTRAC - 1, &
      TAU / 2._dp,TAU,visc)
  CALL output_real_data('vcmac.dat',MAXVCMAC,0._dp,TAU,vcmac)
  CALL output_real_data('testbeta.dat',MAXCHECK,CHECKSTEPS * TAU, &
      CHECKSTEPS * TAU,testbeta)
  CALL output_matrix_data('vcm.dat',VECFORM,MAXCHECK,DIM, &
     CHECKSTEPS * TAU,CHECKSTEPS * TAU,vcm)

CONTAINS

!********************************************************************
  SUBROUTINE init_cell_coords(cvcoord)
!********************************************************************
! This subroutine initializes the 2D coordinates of the collision
! volumes. Note that the coordinates are taken to be the ones in the
! middle of the collision volumes.
!********************************************************************
    
    INTEGER :: &
        ind, &   ! 1D collision volume index
        i, &     ! 2D collision volume index
        j        ! 2D collision volume index
  
    REAL(DP), DIMENSION(1:NWS,1:DIM) :: &
        cvcoord  ! coordinates of collision volumes

!--------------------------------------------------------------------
! To assign the coordinates for the collision volumes, we loop
! over the 2D collision volume indices.
!--------------------------------------------------------------------

    DO i = 1, NWSX
      
      DO j = 1, NWSY

!--------------------------------------------------------------------
! For collision volumes in a lattice l_i x l_j
! with indices i and j we can compute the c.v. indices e.g. 
! as follows: ind = 1 + (j - 1) + l_j * (i - 1).
!--------------------------------------------------------------------

        ind = j + NWSY * (i - 1)        
        cvcoord(ind,X) = LWS * i - LWS / 2._dp
        cvcoord(ind,Y) = LWS * j - LWS / 2._dp

      END DO

    END DO
          
  END SUBROUTINE init_cell_coords


!********************************************************************
  SUBROUTINE init_posit_se(posit,positnpbc)
!********************************************************************
! This subroutine initializes the positions of the solute particles
! to a "random polymer" configuration.
!********************************************************************

    INTEGER :: &
        i         ! particle index 

    REAL(DP) :: &
        ranang    ! random angle

    REAL(DP), DIMENSION(1:NSE - 1) :: &
        ranvect   ! vector of uniformly distributed random numbers

    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        posit     ! particle positions

    REAL(DP), DIMENSION(1:NSE,1:DIM) :: &
        positnpbc ! solute posit. without periodic boundary conditions

!--------------------------------------------------------------------
! The first monomer in the polymer chain in placed, say, in
! the middle of the simulation box.
!--------------------------------------------------------------------

    posit(1,:) = LENG / 2._dp
    positnpbc(1,:) = posit(1,:)

!--------------------------------------------------------------------
! A vector of new uniform random numbers in [0,1] is generated.
!--------------------------------------------------------------------

    CALL unif_ran_array(NSE - 1,1,ranvect)

!--------------------------------------------------------------------
! The rest of the polymer is taken care of. 
!--------------------------------------------------------------------

    DO i = 2, NSE
      
      ranang = (MAXINITANG - MININITANG) * ranvect(i - 1) + MININITANG
      posit(i,X) = posit(i - 1,X) + INITDIST * cos(ranang)
      posit(i,Y) = posit(i - 1,Y) + INITDIST * sin(ranang) 
      positnpbc(i,:) = posit(i,:)

    END DO

!--------------------------------------------------------------------
! Note that the periodic boundaries must be handled afterwards:
! otherwise positnpbc is initialized incorrectly.
!--------------------------------------------------------------------

    DO i = 1, NSE

      posit(i,:) = posit(i,:) - LENG * anint(posit(i,:) / LENG - 0.5_dp)     

    END DO

  END SUBROUTINE init_posit_se


!********************************************************************
  SUBROUTINE init_posit_st(posit)
!********************************************************************
! This subroutine initializes the positions of the solvent particles
! to a random configuration.
!********************************************************************

    INTEGER :: &
        i, &        ! particle index 
        j           ! particle index 

    REAL(DP) :: &
        distijsq, & ! distance between i and j squared
        mindistsq   ! minimum distance between solvent and polymer sq.

    REAL(DP), DIMENSION(1:DIM) :: &
        distij      ! distance between particles i and j

    REAL(DP), DIMENSION(1:NST,1:DIM) :: &
        ranarray    ! array of uniformly distributed random numbers

    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        posit       ! particle positions

!--------------------------------------------------------------------
! An array of new uniform random numbers in [0,1] is generated.
!--------------------------------------------------------------------

    CALL unif_ran_array(NST,DIM,ranarray)

!--------------------------------------------------------------------
! Random values for the coordinates of the solvent particles are
! assigned. Note that a coordinate may assume any value [0,L).
! To make sure that no coordinate is assigned the value L,
! we use periodic boundary conditions. 
!--------------------------------------------------------------------

    DO i = NSE + 1, N

      posit(i,:) = LENG * ranarray(i - NSE,:)
      posit(i,:) = posit(i,:) - &
          LENG * anint(posit(i,:) / LENG - 0.5_dp)
      
!--------------------------------------------------------------------
! Check that the solvent particle is not too close to a solute
! particle. (In that case there would be enormous potential
! gradients.) If the distance is too short, the solvent
! particle is given a new random position. Note that we need not
! check for solvent-solvent overlaps: these do not matter when we
! are dealing with a coarse-grained solvent.
!--------------------------------------------------------------------

      DO

        mindistsq = SIG * SIG
        
        DO j = 1, NSE

          distij = posit(i,:) - posit(j,:)
          distij = distij - LENG * anint(distij / LENG)
          distijsq = sum(distij * distij)
          mindistsq = min(distijsq,mindistsq)
 
        END DO
        
        IF (mindistsq >= SIG * SIG) EXIT

        posit(i,X) = LENG(X) * grnd()
        posit(i,Y) = LENG(Y) * grnd()
        posit(i,:) = posit(i,:) - &
            LENG * anint(posit(i,:) / LENG - 0.5_dp)
       
      END DO

    END DO

  END SUBROUTINE init_posit_st


!********************************************************************
  SUBROUTINE init_veloc(veloc)
!********************************************************************
! This subroutine initializes the velocities of all particles using 
! a Maxwellian velocity distribution. As solvents and solutes
! may have different masses, they are initialized separately from
! different distributions. 
!********************************************************************

    INTEGER :: &
        i             ! particle index 

    REAL(DP) :: &
        meanse, &     ! mean of Gaussian random numbers
        meanst, &     ! mean of Gaussian random numbers
        variancese, & ! variance of Gaussian random numbers
        variancest    ! variance of Gaussian random numbers

    REAL(DP), DIMENSION(1:NSE,1:DIM) :: &
        ranarrayse    ! array of Gaussian random numbers
    
    REAL(DP), DIMENSION(1:NST,1:DIM) :: &
        ranarrayst    ! array of Gaussian random numbers

    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        veloc         ! particle velocities

!--------------------------------------------------------------------
! Two array of new Gaussian random numbers are generated: one
! for the solvwnt and another for the solute.
! The means are zero and the variances are (k_B T) / m.
!--------------------------------------------------------------------
    
    meanse = UMEAN
    meanst = UMEAN
    variancese = 1._dp / (MSE * BETA)
    variancest = 1._dp / (MST * BETA)

    CALL gauss_ran_array(NSE,DIM,meanse,variancese,ranarrayse)
    CALL gauss_ran_array(NST,DIM,meanst,variancest,ranarrayst)

    DO i = 1, NSE
      
      veloc(i,:) = ranarrayse(i,:)
  
    END DO
    
    DO i = NSE + 1, N
      
      veloc(i,:) = ranarrayst(i - NSE,:)
  
    END DO

!--------------------------------------------------------------------
! The drift must be removed and the velocities must be scaled to
! obtain the correct temperature.
!--------------------------------------------------------------------

    CALL remove_drift(veloc)
    CALL scale_veloc(veloc)
    
  END SUBROUTINE init_veloc


!********************************************************************
  SUBROUTINE scale_veloc(veloc)    
!********************************************************************
! This subroutine scales the Maxwell distributed velocities so
! that the desired temperature is obtained. 
!********************************************************************

    INTEGER :: &
        i           ! particle index 

    REAL(DP) :: &
        sumsqse, &  ! sum of solute velocities squared
        sumsqst, &  ! sum of solvent velocities squared
        scal        ! scaling factor for velocities

    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        veloc       ! particle velocities

!--------------------------------------------------------------------
! The scaling factor is computed as sqrt(k_B T N d / (sum(m_i v_i^2))),
! and each velocity must be scaled by this factor.
!--------------------------------------------------------------------

    sumsqse = 0._dp
    sumsqst = 0._dp

    DO i = 1, NSE

      sumsqse = sumsqse + sum(veloc(i,:) * veloc(i,:))

    END DO

    DO i = NSE + 1, N

      sumsqst = sumsqst + sum(veloc(i,:) * veloc(i,:))

    END DO
    
    scal = sqrt((N * DIM) / (BETA * (MSE * sumsqse + MST * sumsqst)))
    
    DO i = 1, N

      veloc(i,:) = scal * veloc(i,:)

    END DO

  END SUBROUTINE scale_veloc


!********************************************************************
  SUBROUTINE remove_drift(veloc)
!********************************************************************
! This subroutine removes any drift from the system.
!********************************************************************

    INTEGER :: &
        i, &        ! particle index 
        alpha       ! coordinate index

    REAL(DP), DIMENSION(DIM) :: &
        avese, &    ! average velocity for solute
        avest       ! average velocity for solvent

    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        veloc       ! particle velocities
    
!--------------------------------------------------------------------
! The drift is removed. The polymer and the solvent are
! treated separately, i.e. the polymer is not moving.
!--------------------------------------------------------------------

    avese = (/ (0._dp, alpha = 1, DIM) /)    
    avest = (/ (0._dp, alpha = 1, DIM) /)    

    DO i = 1, NSE

      avese = avese + veloc(i,:)

    END DO

    DO i = NSE + 1, N
      
      avest = avest + veloc(i,:)

    END DO

    avese = avese / NSE
    avest = avest / NST

    DO i = 1, NSE

      veloc(i,:) = veloc(i,:) - avese

    END DO

    DO i = NSE + 1, N

      veloc(i,:) = veloc(i,:) - avest

    END DO

  END SUBROUTINE remove_drift


!*******************************************************************
  SUBROUTINE init_ac(nac,maxt,tac,account,tacbound,acloop,ac)
!*******************************************************************
! This subroutine initializes the quantities needed for the
! "on the fly" computation of an autocorrelation function.
!*******************************************************************

    INTEGER :: &
        nac, &      ! number of autocorrelation functions
        maxt, &     ! size of autocorrelation function
        tac, &      ! time index for autocorrelation function
        account, &  ! counter for autocorrelation function
        it, &       ! time index
        ic          ! a.c. function index

    INTEGER, DIMENSION(1:maxt + 1) :: &
        tacbound    ! for calculating a.c. f. "on the fly"
  
    INTEGER, DIMENSION(-maxt:maxt) :: &
        acloop      ! for calculating a.c. f. "on the fly"

    REAL(DP), DIMENSION(1:nac,0:maxt - 1) :: &
        ac          ! autocorrelation function
   
    tac = 1
    account = 0
    tacbound = (/ (it, it = 1, maxt), 1 /)
    acloop = (/ (it, it = 0, maxt), (it, it = 1, maxt) /)

    DO ic = 1, nac
      
      ac(ic,:) = (/ (0._dp, it = 1, maxt) /)
      
    END DO

  END SUBROUTINE init_ac


!*******************************************************************
  SUBROUTINE assign_cell_indices(cvind,posit)
!*******************************************************************
! This subroutine assigns collision volume indices 
! for the solvent particles. (For solutes the indices are never
! used.)
!*******************************************************************

    INTEGER :: &
        i        ! particle index

    INTEGER, DIMENSION(1:NST) :: &
        cvind    ! indicates to which coll. vol. a particle belongs

    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        posit    ! particle positions

    DO i = NSE + 1, N

!--------------------------------------------------------------------
! The collision volume index for a given particle is computed.
! For collision volumes in a lattice l_i x l_j
! with indices i and j we can compute the c.v. indices e.g. 
! as follows: ind = 1 + (j - 1) + l_j * (i - 1).
!--------------------------------------------------------------------
      
      cvind(i - NSE) = 1 + int(posit(i,Y) / LWS) + &
          NWSY * int(posit(i,X) / LWS)

    END DO

  END SUBROUTINE assign_cell_indices


!********************************************************************
  SUBROUTINE unif_ran_array(rows,columns,ranarray)
!********************************************************************
! This subroutine produces a 2D array containing pseudorandom 
! real numbers (double). These are uniformly distributed and [0,1].
!********************************************************************

    INTEGER :: &
        rows, &    ! number of rows in the array
        columns, & ! number of columns in the array
        i, &       ! row index
        j          ! column index
    
    REAL(DP), DIMENSION(1:rows,1:columns) :: &
        ranarray   ! array containing random numbers

    DO i = 1, rows

      DO j = 1, columns
      
        ranarray(i,j) = grnd()
      
      END DO

    END DO

END SUBROUTINE unif_ran_array


!********************************************************************
  SUBROUTINE gauss_ran_array(rows,columns,mean,variance,ranarray)
!********************************************************************
! This subroutine produces a 2D array containing pseudorandom 
! real numbers (double). These have a Gaussian distribution
! and a mean and a variance specified by the user. 
!********************************************************************

    INTEGER :: &
        rows, &     ! number of rows in the array
        columns, &  ! number of columns in the array
        last, &     ! last row to be considered in the loop over rows
        i, &        ! row index
        j           ! column index
    
    REAL(DP) :: &
        mean, &     ! mean of Gaussian random numbers
        variance, & ! variance of Gaussian random numbers
        fact1, &    ! common factor in two Gaussian random numbers
        fact2       ! common factor in two Gaussian random numbers

    REAL(DP), DIMENSION(1:rows,1:columns) :: &
        ranarray    ! array containing random numbers

!--------------------------------------------------------------------
! An array of new uniform random numbers in [0,1] is generated.
!--------------------------------------------------------------------

    CALL unif_ran_array(rows,columns,ranarray)
    
!--------------------------------------------------------------------
! The uniform random numbers are converted to Gaussian ones
! with the desired properties. This is done using the Box-Muller
! algorithm decribed in A & T, p. 347. The algorithm
! produces *pairs* of Gaussian distributed random numbers,
! hence the peculiar DO and IF loops etc. (The number of 
! elements might be odd...)
!--------------------------------------------------------------------

    IF (modulo(rows,2) == 0) THEN
      
      last = rows - 1
    
    ELSE
    
      last = rows - 2
    
    END IF

    DO j = 1, columns
    
      DO i = 1, last, 2
    
        fact1 = sqrt(-2._dp * variance * log(ranarray(i,j)))
        fact2 = TWOPI_D * ranarray(i + 1,j)
        ranarray(i,j) = mean + fact1 * cos(fact2)         
        ranarray(i + 1,j) = mean + fact1 * sin(fact2)
      
      END DO

      IF (modulo(rows,2) /= 0) THEN
      
        fact1 = sqrt(-2._dp * variance * log(ranarray(rows,j)))
        fact2 = TWOPI_D * grnd()
        ranarray(rows,j) = mean + fact1 * cos(fact2)
      
      END IF
    
    END DO

END SUBROUTINE gauss_ran_array


!*******************************************************************
  SUBROUTINE collide_st(cvamount,cvaveveloc,cvrotation,cvind,veloc)
!*******************************************************************
! This subroutine handles the collision part of the solvent dynamics. 
!*******************************************************************

    INTEGER :: &
        i, &        ! particle index
        ind         ! collision volume index

    REAL(DP) :: &
        ranang, &   ! random angle
        cra, &      ! cosine of random angle
        sra         ! cosine of random angle

    REAL(DP), DIMENSION(1:NWS) :: &
        ranvect     ! vector containing random numbers
    
    REAL(DP), DIMENSION(1:NWS) :: &
        cvamount    ! amount of particles in collision volume

    REAL(DP), DIMENSION(1:NWS,1:DIM) :: &
        cvaveveloc  ! average velocities of collision volumes

    REAL(DP), DIMENSION(1:NWS,1:DIM,1:DIM) :: &
        cvrotation  ! rotation matrices for collision volumes

    INTEGER, DIMENSION(1:NST) :: &
        cvind       ! indicates to which coll. vol. a particle belongs

    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        veloc       ! particle velocities

!--------------------------------------------------------------------
! Initialize the mean velocities and particle numbers in the
! collision volumes to zero.
!--------------------------------------------------------------------

    DO ind = 1, NWS

      cvaveveloc(ind,:) = (/ (0._dp, ind = 1, DIM) /)
      cvamount(ind) = 0

    END DO

!--------------------------------------------------------------------
! Compute the average velocity component and 
! total amount of particles for each collision volume.
! Notice that only solvent particles participate in collisions.
!--------------------------------------------------------------------

    DO i = NSE + 1, N

!--------------------------------------------------------------------
! The collision volume indices of the particles have already
! been updated. Now they are being used for identifying the
! appropriate collision volume for each particle.
!--------------------------------------------------------------------
      
      ind = cvind(i - NSE)

!--------------------------------------------------------------------
! After identifying the cell, we increment the appropriate
! velocity components and the particle number in that cell.      
!--------------------------------------------------------------------

      cvaveveloc(ind,:) = &
          cvaveveloc(ind,:) + veloc(i,:)

     cvamount(ind) = cvamount(ind) + 1

    END DO

!--------------------------------------------------------------------
! The average velocities must be scaled by the number of particles
! in that volume. Note that the particle number in a box might 
! well be zero. In that case the average velocity is automatically
! zero due to initialization.
!--------------------------------------------------------------------

    DO ind = 1, NWS

      IF (cvamount(ind) /= 0) THEN

        cvaveveloc(ind,:) = cvaveveloc(ind,:) / cvamount(ind)

      END IF

    END DO

    
!--------------------------------------------------------------------
! An array of new uniform random numbers in [0,1] is generated.
!--------------------------------------------------------------------
    
    CALL unif_ran_array(NWS,1,ranvect)

!--------------------------------------------------------------------
! The random rotation matrices for the different cells are computed.
!--------------------------------------------------------------------

    DO ind = 1, NWS

!--------------------------------------------------------------------
! The angle that determines the rotation is chosen randomly.
!--------------------------------------------------------------------

      ranang =  TWOPI_D * ranvect(ind)

!--------------------------------------------------------------------
! Trigonometric functions for the selected angle are "computed".
!--------------------------------------------------------------------

      cra = cos(ranang) 
      sra = sin(ranang) 

!--------------------------------------------------------------------
! The rotation matrix is computed. See e.g. Goldstein.
!--------------------------------------------------------------------

      cvrotation(ind,X,X) = cra
      cvrotation(ind,X,Y) = sra
      cvrotation(ind,Y,X) = -sra
      cvrotation(ind,Y,Y) = cra

    END DO

!--------------------------------------------------------------------
! The velocities of the individual particles are rotated.
!--------------------------------------------------------------------

    DO i = NSE + 1, N      

      ind = cvind(i - NSE)
      veloc(i,:) = cvaveveloc(ind,:) + matmul(cvrotation(ind,:,:), &
          veloc(i,:) - cvaveveloc(ind,:))

    END DO

  END SUBROUTINE collide_st


!*******************************************************************
  SUBROUTINE do_md(point,list,poten,virial,acc,posit,veloc,positnpbc)
!*******************************************************************
! This subroutine takes care of the MD part of the dynamics.
! The "Velocity Verlet" algorithm is used as an integrator.
!*******************************************************************

    LOGICAL :: &
        update      ! if true the list is updated

    INTEGER :: &
        i           ! particle index

    INTEGER, DIMENSION(1:NSE + 1) :: &
        point       ! index to the Verlet neighbor list
    
    INTEGER, DIMENSION(1:MAXNAB) :: &
        list        ! Verlet neighbour list

    REAL(DP) :: &
        poten, &    ! total potential energy of the system
        virial      ! virial for the system

    REAL(DP), DIMENSION(1:DIM) :: &
        displ       ! displacement of particle during time step

    REAL(DP), DIMENSION(1:NSE,1:DIM) :: &
        positnpbc   ! solute posit. without periodic boundary conditions
    
    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        acc, &      ! particle accelerations
        posit, &    ! particle positions
        veloc, &    ! particle velocities
        oldacc      ! accelerations before the update of positions
    
!--------------------------------------------------------------------
! The accelerations computed at the previous step are stored.
!--------------------------------------------------------------------

    oldacc =  acc

!--------------------------------------------------------------------
! The positions of the particles are updated. Notice that the
! solute positions without the periodic boundary conditions must
! be treated separately.
!--------------------------------------------------------------------

    DO i = 1, NSE

      displ = DT * veloc(i,:) + DT * DT * acc(i,:) / 2._dp
      posit(i,:) = posit(i,:) + displ
      positnpbc(i,:) = positnpbc(i,:) + displ

    END DO

    DO i = NSE + 1, N

      posit(i,:) = posit(i,:) + &
          DT * veloc(i,:) + DT * DT * acc(i,:) / 2._dp

    END DO
  
!--------------------------------------------------------------------
! The accelerations are computed here.  First we need to check whether
! or not the Verlet list should be updated.
!--------------------------------------------------------------------

    CALL check_list(update,oldposit,posit)
    CALL comp_acc(update,point,list,poten,virial,acc,posit)

!--------------------------------------------------------------------
! The velocities of all particles are updated.
!--------------------------------------------------------------------

    DO i = 1, N
      
      veloc(i,:) = veloc(i,:) + &
          DT * (oldacc(i,:) + acc(i,:)) / 2._dp

    END DO
    
  END SUBROUTINE do_md


!*******************************************************************
  SUBROUTINE fold_back(posit)
!*******************************************************************
! This subroutine folds back the particles to the original
! simulation box.
!*******************************************************************

    INTEGER :: &
        i           ! particle index

    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        posit       ! particle positions
    
    DO i = 1, N

      posit(i,:) = posit(i,:) - LENG * anint(posit(i,:) / LENG - 0.5_dp)

    END DO
    
  END SUBROUTINE fold_back


!*******************************************************************
  SUBROUTINE comp_acc(update,point,list,poten,virial,acc,posit)
!*******************************************************************
! This subroutine computes the accelerations.
!*******************************************************************

    LOGICAL :: &
        update      ! if true the list is updated
    
    INTEGER :: &
        i, &         ! particle index
        j, &         ! particle index
        alpha, &     ! coordinate index
        nlist, &     ! list index
        jbeg, &      ! first neighbor of i
        jend, &      ! last neighbor of i
        jnab         ! list index

    INTEGER, DIMENSION(1:NSE + 1) :: &
        point        ! index to the Verlet neighbor list
    
    INTEGER, DIMENSION(1:MAXNAB) :: &
        list         ! Verlet neighbour list

    REAL(DP) :: &
        poten, &     ! total potential energy of the system
        virial, &    ! virial for the system
        distijsq, &  ! distance between i and j squared
        distijmod, & ! modulus of distance between i and j        
        sr2, &       ! (\sigma / r_ij )^2, useful in LJ-pot.
        sr6, &       ! (\sigma / r_ij )^6, useful in LJ-pot.
        potenij, &   ! pair potential for i and j
        virialij     ! contribution to virial from for i and j

    REAL(DP), DIMENSION(1:DIM) :: &
        positi, &    ! position of the ith particle, temp. variable
        acci, &      ! acceleration of the ith particle, temp. variable
        distij, &    ! distance between particles i and j
        accij        ! force on j due to i (not yet acceleration)
   
    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        acc, &       ! particle accelerations
        posit        ! particle positions
 
!--------------------------------------------------------------------
! The accelerations, potential energy, and virial are initialized 
! to zero.
!--------------------------------------------------------------------    

    DO i = 1, N

      acc(i,:) = (/ (0._dp, alpha = 1, DIM) /)

    END DO
    
    poten = 0._dp
    virial = 0._dp

!--------------------------------------------------------------------
! First, the truncated Lennard-Jones interactions are considered. 
! For this we need a Verlet neighbor list. 
! If the list is to be updated, we update it and compute the forces. 
!--------------------------------------------------------------------      

    IF (update) THEN

!--------------------------------------------------------------------
! The current configuration is saved.
!--------------------------------------------------------------------
      
      oldposit = posit

      nlist = 0

      DO i = 1, NSE

        point(i) = nlist + 1

!--------------------------------------------------------------------
! The use of temporary variables saves time, as positions and 
! accelerations need no longer be fetched many times
! from large arrays of structs.
!--------------------------------------------------------------------

        positi = posit(i,:)
        acci = acc(i,:) 

        DO j = i + 1, N
          
          distij = positi - posit(j,:)
          distij = distij - LENG * anint(distij / LENG)
          distijsq = sum(distij * distij)
          
          IF (distijsq < RLIST * RLIST) THEN
            
            nlist = nlist + 1
            list(nlist) = j
            
!--------------------------------------------------------------------
! Check if MAXNAB is large enough. (This can be commented out later.)
!--------------------------------------------------------------------

            IF (nlist == MAXNAB) THEN

              stop 'Verlet list too small.'

            END IF
            
            IF (distijsq < RCUT * RCUT) THEN

              sr2 = SIG * SIG / distijsq
              sr6 = sr2 * sr2 * sr2
              potenij = sr6 * (sr6 - 1._dp) + 0.25_dp
              virialij = sr6 * (sr6 - 0.5_dp)
              poten = poten + 4._dp * EPS * potenij
              virial = virial + 48._dp * EPS * virialij / DIM
              accij = 48._dp * EPS * virialij * distij / distijsq
              acci = acci + accij
              acc(j,:) = acc(j,:) - accij
              
            END IF
             
          END IF

        END DO

        acc(i,:) = acci

        END DO

        point(NSE + 1) = nlist + 1


!--------------------------------------------------------------------
! If the list need not be updated, we just use the old list to
! find the neighbors and compute the Lennard-Jones type forces.
!--------------------------------------------------------------------

      ELSE

        DO i = 1, NSE

          jbeg = point(i)
          jend = point(i + 1) - 1

!--------------------------------------------------------------------
! First we must check that the particle really has neighbors.
!--------------------------------------------------------------------

          IF (jbeg <= jend) THEN

            positi = posit(i,:)
            acci = acc(i,:)

            DO jnab = jbeg, jend

              j = list(jnab)

              distij = positi - posit(j,:)
              distij = distij - LENG * anint(distij / LENG)
              distijsq = sum(distij * distij)

              IF (distijsq < RCUT * RCUT) THEN

              sr2 = SIG * SIG / distijsq
              sr6 = sr2 * sr2 * sr2
              potenij = sr6 * (sr6 - 1._dp) + 0.25_dp
              virialij = sr6 * (sr6 - 0.5_dp)
              poten = poten + 4._dp * EPS * potenij
              virial = virial + 48._dp * EPS * virialij / DIM
              accij = 48._dp * EPS * virialij * distij / distijsq
              acci = acci + accij
              acc(j,:) = acc(j,:) - accij
              
            END IF

          END DO

          acc(i,:) = acci

        END IF

      END DO

    END IF

!--------------------------------------------------------------------
! The spring potentials for nearest-neighbor solutes are considered.
!--------------------------------------------------------------------

    DO i = 1, NSE - 1
      
      distij = posit(i,:) - posit(i + 1,:)
      distij = distij - LENG * anint(distij / LENG)
      distijsq = sum(distij * distij)
      distijmod = sqrt(distijsq)

      potenij = (distijmod - LB) * (distijmod - LB) 
      virialij = (distijsq - LB * distijmod)
      poten = poten + KB * potenij / 2._dp
      virial = virial - KB * virialij / DIM 
      accij = -KB * virialij * distij / distijsq
      acc(i,:) = acc(i,:) + accij
      acc(i + 1,:) = acc(i + 1,:) - accij

    END DO

!--------------------------------------------------------------------
! The accelerations are scaled.
!--------------------------------------------------------------------

    DO i = 1, NSE
      
      acc(i,:) = acc(i,:) / MSE

    END DO

    DO i = NSE + 1, N
      
      acc(i,:) = acc(i,:) / MST

    END DO

  END SUBROUTINE comp_acc


!*******************************************************************
  SUBROUTINE check_list(update,oldposit,posit)
!*******************************************************************
! This subroutine indicates whether the Verlet list needs to be
! reconstructed.
!*******************************************************************

    LOGICAL :: &
        update      ! if true the list is updated

    INTEGER :: &
        i           ! particle index
   
    REAL(DP) :: &
        maxdisp     ! the largest displacement

    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        oldposit, & ! particle positions
        posit       ! particle positions

!--------------------------------------------------------------------
! Initialize the largest displacement to zero.
!--------------------------------------------------------------------    
   
    maxdisp = 0._dp

!--------------------------------------------------------------------
! The largest displacement since the last update of the Verlet
! list is computed. 
! Test if this can be done in a more elegant manner:
! maxdisp(X,Y,Z,maxdisp).
!--------------------------------------------------------------------    

    DO i = 1, N

      maxdisp = max(abs(posit(i,X) - oldposit(i,X)),maxdisp)
      maxdisp = max(abs(posit(i,Y) - oldposit(i,Y)),maxdisp)

    END DO

!--------------------------------------------------------------------
! A conservative test for the skin crossing is performed.
!--------------------------------------------------------------------

    maxdisp = 2._dp * sqrt(DIM * maxdisp * maxdisp)

    update = (maxdisp > (RLIST - RCUT))

  END SUBROUTINE check_list


!*******************************************************************
  SUBROUTINE comp_deltaxii(cvcoord,cvind,deltaxii,posit)
!*******************************************************************
! This subroutine computes the displacement in the y direction 
! from the previous collision volume to the current one.
! The displacement is only computed for the solvent.
!*******************************************************************

    INTEGER :: &
        i        ! particle index

    REAL(DP), DIMENSION(1:NWS,1:DIM) :: &
        cvcoord  ! coordinates of collision volumes

    INTEGER, DIMENSION(1:NST) :: &
        cvind    ! indicates to which coll. vol. a particle belongs
  
    REAL(DP), DIMENSION(1:NST) :: &
        deltaxii ! \xi_iy(\tau + 1) - \xi_iy(\tau)

    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        posit    ! particle positions

    DO i = NSE + 1, N

      deltaxii(i - NSE) =  LWS * (1._dp + floor(posit(i,Y) / LWS)) - &
          LWS / 2._dp - cvcoord(cvind(i - NSE),Y)

    END DO

  END SUBROUTINE comp_deltaxii


!*******************************************************************
  SUBROUTINE comp_stress(str,deltaxii,veloc)
!*******************************************************************
! This subroutine is used for computing the elements
! needed for the stress-stress autocorrelation function
! defined in Malevanets et al. Only the solvent particles
! are included in this computation.
!*******************************************************************

    INTEGER :: &
        i           ! particle index 

    REAL(DP) :: &
        str         ! element for the stress-stress autocorrelation (M.)

    REAL(DP), DIMENSION(1:NST) :: &
        deltaxii    ! \xi_iy(\tau + 1) - \xi_iy(\tau)

    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        veloc       ! particle velocities

!--------------------------------------------------------------------    
! For the definition of the stress-stress autocorrelation
! function see Malevanets et al., J. Chem. Phys. 112, 7260 (2000).
!--------------------------------------------------------------------    
 
    str = 0._dp

    DO i = NSE + 1, N

      str = str + veloc(i,X) * deltaxii(i - NSE)

    END DO

  END SUBROUTINE comp_stress


!*******************************************************************
  SUBROUTINE comp_ac(nac,maxt,datadim,tac,account,tacbound, &
      acloop,ac,acstore)
!*******************************************************************
! This subroutine handles the "on the fly" computation of an
! autocorrelation function. This subroutine is redundant: 
! the same thing could be achieved by calling comp_cc with 
! suitable arguments. However, not having to transfer the same 
! array twice might save time.
!*******************************************************************
   
    INTEGER :: &
        nac, &      ! number of autocorrelation functions
        maxt, &     ! size of autocorrelation function
        datadim, &  ! dimension of the stored data elements
        tac, &      ! time index for autocorrelation function
        account, &  ! counter for autocorrelation function
        tdiff, &    ! difference between two times
        hist, &     ! time index for computing the a.c. f.
        ic          ! a.c. function index

    INTEGER, DIMENSION(1:maxt + 1) :: &
        tacbound    ! for calculating a.c. f. "on the fly"
  
    INTEGER, DIMENSION(-maxt:maxt) :: &
        acloop      ! for calculating a.c. f. "on the fly"    

    REAL(DP), DIMENSION(1:nac,0:maxt - 1) :: &
        ac          ! autocorrelation function
    
    REAL(DP), DIMENSION(1:nac,1:maxt,1:datadim) :: &
        acstore     ! stores maxt quantities for computation of a.c. f.

!--------------------------------------------------------------------
! The counter is updated.
!--------------------------------------------------------------------    

    account = account + 1

!--------------------------------------------------------------------
! Until the storing array becomes full for the first time,
! we do nothing but storing. Then we may start incrementing
! the autocorrelation function. 
!--------------------------------------------------------------------

    IF (account >= maxt) THEN

      DO tdiff = 0, maxt - 1

        hist = acloop(tac - tdiff)
        
        DO ic = 1, nac

          ac(ic,tdiff) = ac(ic,tdiff) + &
              sum(acstore(ic,tac,:) * acstore(ic,hist,:))

        END DO

      END DO
 
    END IF

!--------------------------------------------------------------------
! The "storing index" is updated.
!--------------------------------------------------------------------

    tac = tacbound(tac + 1)

  END SUBROUTINE comp_ac


!*******************************************************************
  SUBROUTINE fin_ac(nac,maxt,account,ac)
!*******************************************************************
! This subroutine completes the "on the fly" computation of the
! autocorrelation function. (Can also be used for cross correlation
! functions and mean squares.)
!*******************************************************************
   
    INTEGER :: &
        nac, &     ! number of autocorrelation functions
        maxt, &    ! size of autocorrelation function
        account    ! counter for autocorrelation function

    REAL(DP), DIMENSION(1:nac,0:maxt - 1) :: &
        ac         ! autocorrelation function

    IF (account > maxt) THEN

      ac = ac / real(account - maxt + 1,kind = DP)
      
    END IF

  END SUBROUTINE fin_ac


!*******************************************************************
  SUBROUTINE comp_visc(strac,visc)
!*******************************************************************
! This subroutine is used for computing the viscosity from the
! stress-stress autocorrelation function.
!*******************************************************************

    INTEGER :: &
        it          ! time index 

    REAL(DP) :: &
        scal        ! scaling factor for viscosity

    REAL(DP), DIMENSION(0:MAXSTRAC - 1) :: &
        strac       ! stress-stress autocorrelation function (M.)

    REAL(DP), DIMENSION(0:MAXSTRAC - 2) :: &
        visc        ! viscosity (Malevanets)
        
!--------------------------------------------------------------------    
! For computation of viscosity, see Malevanets et al., 
! J. Chem. Phys. 112, 7260 (2000). Basically the stress-stress
! autocorrelation function is integrated over time and scaled.
! The trapezoid approximation is used.
!--------------------------------------------------------------------    

!--------------------------------------------------------------------    
! The scaling factor is computed. The time step and mass squared 
! must be included---in Malevanets et al. these are always unity.
!--------------------------------------------------------------------    

    scal = MST * MST * RHO * BETA / (2._dp * NST * TAU)  

!--------------------------------------------------------------------    
! The viscosity is initialized.
!--------------------------------------------------------------------    
    
    visc = (/ (0._dp, it = 1, MAXSTRAC - 1) /)
    visc(0) = scal * (strac(0) + strac(1))

!--------------------------------------------------------------------    
! The viscosity is computed and scaled.
!--------------------------------------------------------------------    

    DO it = 1, MAXSTRAC - 2

      visc(it) = visc(it - 1) + scal * (strac(it) + strac(it + 1)) 
      
    END DO

  END SUBROUTINE comp_visc


!*******************************************************************
  SUBROUTINE check_temp(testbeta,veloc)
!*******************************************************************
! This subroutine is used for checking the temperature of the 
! system.
!*******************************************************************

    INTEGER :: &
        i           ! particle index 
    
    REAL(DP) :: &
        testbeta, & ! computed value of 1 / (k_B T)
        sumsqse, &  ! sum of solute velocities squared
        sumsqst     ! sum of solvent velocities squared

    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        veloc       ! particle velocities

!--------------------------------------------------------------------    
! 1 / (k_B T) is calculated from the kinetic energy of the 
! particles as follows: N d / sum(m_i v_i^2).
!--------------------------------------------------------------------    

    sumsqse = 0._dp
    sumsqst = 0._dp

    DO i = 1, NSE

      sumsqse = sumsqse + sum(veloc(i,:) * veloc(i,:))

    END DO

    DO i = NSE + 1, N

      sumsqst = sumsqst + sum(veloc(i,:) * veloc(i,:))

    END DO

    testbeta = N * DIM / (MSE * sumsqse + MST * sumsqst)

  END SUBROUTINE check_temp


!*******************************************************************
  SUBROUTINE comp_centmass(centmass,positnpbc)
!*******************************************************************
! This subroutine computes the center of mass of the polymer.
!*******************************************************************

    INTEGER :: &
        i, &        ! particle index
        alpha       ! coordinate index

    REAL(DP), DIMENSION(1:DIM) :: &
        centmass    ! center of mass of the polymer

    REAL(DP), DIMENSION(1:NSE,1:DIM) :: &
        positnpbc   ! solute posit. without periodic boundary conditions
    
    centmass = (/ (0._dp, alpha = 1, DIM) /)

    DO i = 1, NSE

      centmass = centmass + positnpbc(i,:)

    END DO

    centmass = centmass / NSE

  END SUBROUTINE comp_centmass


!*******************************************************************
  SUBROUTINE comp_radgyr(radgyrsq,centmass,positnpbc)
!*******************************************************************
! This subroutine computes the radius of gyration of the polymer.
!*******************************************************************

    INTEGER :: &
        i           ! particle index

    REAL(DP) :: &
        radgyrsq    ! radius of gyration of the polymer squared

    REAL(DP), DIMENSION(1:DIM) :: &
        centmass, & ! center of mass of the polymer
        dist        ! distance between particle i and center of mass

    REAL(DP), DIMENSION(1:NSE,1:DIM) :: &
        positnpbc   ! solute posit. without periodic boundary conditions

    radgyrsq = 0._dp

    DO i = 1, NSE
      
      dist = positnpbc(i,:) - centmass
      radgyrsq = radgyrsq + sum(dist * dist)

    END DO

    radgyrsq = radgyrsq / NSE

  END SUBROUTINE comp_radgyr


!*******************************************************************
  SUBROUTINE comp_endtoend(endtoendsq,positnpbc)
!*******************************************************************
! This subroutine computes the end-to-end distance of the polymer.
!*******************************************************************

    INTEGER :: &
        i           ! particle index

    REAL(DP) :: &
        endtoendsq  ! end-to-end distance of the polymer

    REAL(DP), DIMENSION(1:DIM) :: &
        dist        ! distance between particles 1 and NSE

    REAL(DP), DIMENSION(1:NSE,1:DIM) :: &
        positnpbc   ! solute posit. without periodic boundary conditions

    dist = positnpbc(NSE,:) - positnpbc(1,:)
    endtoendsq = sum(dist * dist)

  END SUBROUTINE comp_endtoend


!*******************************************************************
  SUBROUTINE comp_cmveloc(cmveloc,veloc)
!*******************************************************************
! This subroutine computes the center of mass velocity of the
! polymer.
!*******************************************************************

    INTEGER :: &
        i, &        ! particle index
        alpha       ! coordinate index

    REAL(DP), DIMENSION(1:DIM) :: &
        cmveloc     ! center of mass velocity of the polymer

    REAL(DP), DIMENSION(1:N,1:DIM) :: &
        veloc       ! particle velocities
    
    cmveloc = (/ (0._dp, alpha = 1, DIM) /)

    DO i = 1, NSE

      cmveloc = cmveloc + veloc(i,:)

    END DO

    cmveloc = cmveloc / NSE

  END SUBROUTINE comp_cmveloc


!*******************************************************************
  SUBROUTINE output_configuration(filename,form,positnpbc)
!*******************************************************************
! This subroutine outputs the coordinates of the monomers.
!*******************************************************************

    INTEGER :: &
        i, &         ! particle index
        alpha        ! coordinate index
    
    CHARACTER(LEN = *) :: &
        form         ! file format

    CHARACTER(LEN = MAXSTR) :: &
        filename     ! the name of the output file 

    REAL(DP), DIMENSION(1:NSE,1:DIM) :: &
        positnpbc    ! solute posit. without periodic boundary conditions
  
!--------------------------------------------------------------------
! The output file is opened.
!--------------------------------------------------------------------
    
    open(unit = 1,status = 'new',file = filename, &
        form = 'formatted', access = 'sequential')

!--------------------------------------------------------------------
! The output data is written.
!--------------------------------------------------------------------    

    write(1,*) '# ', NSE
    write(1,*) '# Comments'

    DO i = 1, NSE
      
      write(1,form) positnpbc(i,:)

    END DO

!--------------------------------------------------------------------
! The output file is closed.
!--------------------------------------------------------------------        

    close(1)

  END SUBROUTINE output_configuration


!*******************************************************************
  SUBROUTINE output_general(filename)
!*******************************************************************
! This subroutine outputs general characteristics for the run:
! the number of particles, number of steps etc.
!*******************************************************************
    
    CHARACTER(LEN = MAXSTR) :: &
        filename    ! the name of the output file 
    
!--------------------------------------------------------------------
! The output file is opened.
!--------------------------------------------------------------------
    
    open(unit = 1,status = 'new', file = filename, &
        form = 'formatted', access = 'sequential')

!--------------------------------------------------------------------
! The output data is written.
!--------------------------------------------------------------------    


    write(1,*) 'no. of solutes: ', NSE
    write(1,*) 'no. of solvents: ', NST
    write(1,*) 'no. of coll. vol. in x dir.: ', NWSX
    write(1,*) 'no. of coll. vol. in y dir.: ', NWSY
    write(1,*) 'dimension: ', DIM
    write(1,*) 'random seed: ', RSEED
    write(1,*) 'sigma: ', SIG
    write(1,*) 'epsilon: ', EPS
    write(1,*) 'solvent mass: ', MST
    write(1,*) 'LJ time: ', TAULJ
    write(1,*) 'beta: ', BETA
    write(1,*) 'MD time step: ', DT
    write(1,*) 'collision time step: ', TAU
    write(1,*) 'solute mass: ', MSE
    write(1,*) 'mean velocity of system: ', UMEAN
    write(1,*) 'cutoff radius: ', RCUT
    write(1,*) 'linear size, X: ', LENG(X)
    write(1,*) 'linear size, Y: ', LENG(Y)
    write(1,*) 'tempering steps: ', TEMPSTEPS
    write(1,*) 'equilibration steps: ', EQSTEPS
    write(1,*) 'simulation steps: ', COLSTEPS
    write(1,*) 'collision volume size: ', LWS
    write(1,*) 'density: ', RHO

!--------------------------------------------------------------------
! The output file is closed.
!--------------------------------------------------------------------        

    close(1)

  END SUBROUTINE output_general

!*******************************************************************
  SUBROUTINE output_real_data(filename,elements,initval,incr,data)
!*******************************************************************
! This subroutine outputs a vector containing double precision
! real data to a file with a desired name. The data is output as
! a function of time/distance/etc.
!*******************************************************************

    INTEGER :: &
        elements, & ! number of elements in a certain data
        id          ! data index

    REAL(DP) :: &
        initval, &  ! value of first data point
        incr        ! increment between data points

    REAL(DP), DIMENSION(1:elements) :: &
        data        ! the data to be output

    CHARACTER(LEN = MAXSTR) :: &
        filename    ! the name of the output file 
    
!--------------------------------------------------------------------
! The output file is opened.
!--------------------------------------------------------------------
    
    open(unit = 1,status = 'new',file = filename, &
        form = 'formatted', access = 'sequential')

!--------------------------------------------------------------------
! The output data is written.
!--------------------------------------------------------------------    

    DO id = 0, elements - 1
      
      write(1,*) initval + id * incr, data(id + 1)

    END DO

!--------------------------------------------------------------------
! The output file is closed.
!--------------------------------------------------------------------        

    close(1)

  END SUBROUTINE output_real_data


!*******************************************************************
  SUBROUTINE output_matrix_data(filename,form,rows,columns, &
      initval,incr,data)
!*******************************************************************
! This subroutine outputs a vector containing double precision
! real data to a file with a desired name. The data is output as
! a function of time/distance/etc.
!*******************************************************************

    INTEGER :: &
        rows, &     ! number of rows in a certain data
        columns, &  ! number of columns in a certain data
        ir, &       ! row index
        ic          ! column index

    CHARACTER(LEN = *) :: &
        form        ! file format

    REAL(DP) :: &
        initval, &  ! value of first data point
        incr        ! increment between data points

    REAL(DP), DIMENSION(1:rows,1:columns) :: &
        data        ! the data to be output

    CHARACTER(LEN = MAXSTR) :: &
        filename    ! the name of the output file 
    
!--------------------------------------------------------------------
! The output file is opened.
!--------------------------------------------------------------------
    
    open(unit = 1,status = 'new',file = filename, &
        form = 'formatted', access = 'sequential')

!--------------------------------------------------------------------
! The output data is written.
!--------------------------------------------------------------------    

    DO ir = 0, rows - 1
      
      write(1,form) initval + ir * incr, data(ir + 1,:)

    END DO

!--------------------------------------------------------------------
! The output file is closed.
!--------------------------------------------------------------------        

    close(1)

  END SUBROUTINE output_matrix_data

END PROGRAM mesosim