author: Pawel Jan Zuk
copyright (c) 2021 Pawel Jan Zuk
under GPL v3

***************************
    INTRODUCTION:
***************************


This file contains the open source FORTRAN implementation of the Generalized Rotne-Prager-Yamakawa hydrodynamic interactions in Lees-Edwards periodic boundary conditions with GRPY particle overlaps and some examples how to use them in case of shear flow.

The code calculates 11n x 11n hydrodynamic matrices for n spherical particles with varying radii that can overlap.

In case you use this code, compiled program or code elements, please cite the scientific papers (we understand that it is a lot but it took that much work):
1) The GRPY method
   Wajnryb, E., Mizerski, K.A., Zuk, P.J. and Szymczak, P., 2013. 
   Generalization of the Rotne–Prager–Yamakawa mobility and shear disturbance tensors. 
   Journal of Fluid Mechanics, 731.
2) The polidispersity (different radii)
   Zuk, P.J., Wajnryb, E., Mizerski, K.A. and Szymczak, P., 2014.
   Rotne–Prager–Yamakawa approximation for different-sized particles in application to macromolecular bead models.
   Journal of Fluid Mechanics, 741.
3) Stresslet components in case of calculating e.g. intrinsic viscosity or normal stresses
   Zuk, P.J., Cichocki, B. and Szymczak, P., 2017.
   Intrinsic viscosity of macromolecules within the generalized Rotne-Prager-Yamakawa approximation.
   Journal of Fluid Mechanics, 822.
4) GRPY in Lees-Edwards Periodic Boundary Conditions
   Mizerski, K.A., Wajnryb, E., Zuk, P.J. and Szymczak, P., 2014.
   The Rotne-Prager-Yamakawa approximation for periodic systems in a shear flow.
   The Journal of chemical physics, 140(18), p.184103.
   and
5) this contribution
   Bogdan Cichocki, Piotr Szymczak, and Paweł J. Żuk, 2021.
   Generalized Rotne–Prager–Yamakawa approximation for Brownian dynamics in shear flow in bounded, unbounded, and periodic domains.
   The Journal of Chemical Physics 154, 124905  


***************************
     INSTALLATION:
***************************


For unix systems  with a gfortran compiler a simple Makefile is prepared. It is sufficient to enter the main folder GRPerY and type
make
The executable file will appear in GRPerY/bin folder

For other operating systems and compilers please modify the Makefile in an appropriate fashion.

This code contains all necessary elements to compile and run. A FORTRAN compiler is sufficient. It uses no external libraries, which however can be used to speed up the computations.


***************************
    USAGE INSTRUCTIONS:
***************************

The hydrodynamic matrices are calculated by procedures in 
GRPerY/src/HYDRO/per_GRPY_polyd.f
file. They can be referenced and used as it is shown in a simplistic program for running a simulation, which we included and built around the library.

The program main file is 
GRPerY/src/HYDRO/main.f 
It reads the input and contains main computational loop that invokes time integrating steps and writing the output. There are two Heun type integrating schemes: first without thermal motion
GRPerY/src/STEPPER/stepper_noBrown.f
and second with Brownian motion
GRPerY/src/STEPPER/stepper_Brown.f
In those steppers the procedure calculating the hydrodynamic mobility matrix
GRPERY_MOB(APP,APQ,AQQ,CONF,RADII,NN,LATTICE,EWS_ERR)
is invoked and used for calculating displacements for time integration.

The procedure

GRPERY_MOB(APP,APQ,AQQ,CONF,RADII,NN,LATTICE,EWS_ERR)

has the following arguments:

APP – returns 6NN x 6NN component mobility matrix for transaltional and rotational coupling.
      It is organized with NN x NN blocks that are 6 x 6.
      Each of the blocks contains
        mu_tt(3x3) mu_tr(3x3)
        mu_rt(3x3) mu_rr(3x3) 
      mobility matrix containing interactions between particles i and j. The i and j also denote block numbering.
      For example block 1:6,7:12 has couplings between particles with i=1 and j=2.
APQ – returns 6NN x 5NN component mobility matrix for transaltional and rotational motion coupling with dipolar degrees of freedom.
      It is organized with NN x NN blocks that are 6 x 5. 
      Each of the blocks contains
        mu_td(3x5) 
        mu_rd(3x5) 
      mobility matrix containing interactions between particles i and j. The i and j also denote block numbering. 
      For example block 1:6,6:10 has couplings between particles with i=1 and j=2.
      The dipolar degrees of freedom are presented in the form of 5 linearly independent components (traceless and symmetric). 
      These components are defined in the ./GRPerY/src/TENSORS/tensors.f file inside INIT_Y2() procedure. 
      For example the simple shear component of the U=(gamma z,0,0) flow is given by the 4th cartesian tensor in Y2 matrix. 
      Y2(4,1,3)=-1.D0/SQRT(2.D0)
      Y2(4,3,1)=-1.D0/SQRT(2.D0)
      where first index is coordinate in Y2 matrix and two following indices represent cartesian components.
      A good description of projection tensors in present on page 163 in
           Ekiel-Jeżewska, M.L., Wajnryb, E., Feuillebois, F. and Sellier, A., 2009.
           Precise multipole method for calculating hydrodynamic interactions between spherical particles in the Stokes flow.
           Theoretical methods for micro scale viscous flows, pp.127-172. 
AQQ – returns 5NNx5NN component mobility matrix for coupling of dipolar degrees of freedom.
      It is organized with NN x NN blocks that are 5 x 5.
      Each of the blocks contains
      mu_dd(5x5)  
      mobility matrix containing interactions between particles i and j.  The i and j also denote block numbering. For example block 1:5,6:10 has couplings between particles with i=1 and j=2
CONF – takes 3 x NN matrix with particle positions. Particles do not have to be in the same periodic cell because anytime a distance is calculated it is calculated between the closest neighbors.
RADII – takes NN array of particle radii
NN – takes number of particles
LATTICE – takes 3 x 3 matrix having lattice vectors.
      First index numbers the vector and second, the coordinates of the vector.
      For example lattice of size x=10,y=20,z=30 will have nonzero values
      LATTICE(1,1)=10
      LATTICE(2,2)=20
      LATTICE(3,3)=20
EWS_ERR – takes error for controlling the extent of Ewald summation.
      For detailed information please check eq. 23 and 24 in 
      Mizerski, K.A., Wajnryb, E., Zuk, P.J. and Szymczak, P., 2014.
      The Rotne-Prager-Yamakawa approximation for periodic systems in a shear flow.
      The Journal of chemical physics, 140(18), p.184103.


***************************
        UNITS:
***************************


calculating the mobility matrices we assume that 
\pi \eta = 1
thus the translational mobility coefficient will return
1 / (6 \pi \eta a) = 1 / (6 a)
where a is particle radius.

Other units are specified by the program using the library. 
In the exemplary implementation we have energy unit k_B T. 
Such choice specifies the relation between length scale and time scale: a sphere with radius 0.5 will have a MSD(t) = 2 t.

Please read the source code for more details.


***************************
       EXAMPLES:
***************************


Two examples are present for Brownian and nonBrownian simulations of spheres interacting with Lennard-Jones potential.
The examples are in file GRPerY/examples folder. 
Each of them contains files:
runSimulation – a file containing command to run the example
control_file.dat – a file containing simulation parameters
initial_config.dat -  a file containing initial configuration, size of the box and initial time from which simulation starts, determines lattice skew in the initial time step.


***************************
      REMARKS:
***************************


The implementation uses LAPACK routines for Cholesky decomposition (CHOLESKY subroutine) and matrix inversions (MATREV subroutine) in GRPerY/src/TENSOR/matrices.f. In this implementation all necessary routines are in GRPerY/src/LAPACK/lapblas_double_excerpts.f file.
These are most computational demanding elements of the code. In case of larger computations it is advisable to exchange those implementations to their parallel versions or substitutes.