! *******************************************************************************
! *   Copyright (C) 2007 - 2015 by J.A.Purton                                   *
! *   john.purton[@]stfc.ac.uk                                                  *
! *                                                                             *
! *   Contributors:                                                             *
! *   -------------                                                             *
! *   A.V.Brukhno (C) 2015-2018                                                 *
! *   andrey.brukhno[@]stfc.ac.uk abrukhno[@]gmail.com                          *
! *   - overall optimizations (general MC scheme, energy calculus and MC steps) *
! *   - Free Energy Difference (FED) & order parameters (fed_*_module.f90)      *
! *   - Replica Exchange (RE) & RE-FED extension (rep_exchange_module.f90)      *
! *   - VdW and general potential forms (vdw_*_module.f90 etc)                  *
! *   - User-defined analytical short-ranged VdW forms (vdw_*.f90 & user_*.inc) *
! *   - planar pore 'slit' & 'slab' constraints (slit_module.f90)               *
! *   - planar pore external potentials (external_potential_module.f90)         *
! *   - mean-field approximation (MFA) for long-range electrostatics in slit    *
! *   - refactor of exclusion lists and neighbour lists, incl. optimizations    *
! *   - refactor of atom and molecule selection to avoid 'empty' MC attempts    *
! *   - Cavity-Bias GCMC: optimised sampling on a 3D grid or mesh               *
! *                                                                             *
! *   T.L.Underwood (C) 2015-2020                                               *
! *   t.l.Underwood[@]bath.ac.uk                                                *
! *   - Lattice/Phase-Switch  MC methodology (psmc_*.f90)                       *
! *   - transition matrix (TM) FED MC methodology                               *
! *   - algorithm optimizations (gcmc_*_module.90, MC steps etc)                *
! *   - random number generator debugging & optimization                        *
! *   - FED-GCMC, the counterpart of FED-NPT, samping is done in N instead of V *
! *   - Gibbs Ensemble revision, debugging and optimization                     *
! *   - orientation-dependent potentials (Gay-Berne) (spin_module.f90)          *
! *   - long-range VdW correction                                               *
! *   - general bug fixes                                                       *
! *                                                                             *
! *   R.J.Grant (C) 2016-17                                                     *
! *   r.j.grant[@]bath.ac.uk                                                    *
! *   - extra potential forms (VDW, vdw_direct_module.f90)                      *
! *   - Ewald sums optimization (field.f90, mc_moves.f90, lattice_module.f90)   *
! *   - initial neighbour-list optimization (field.f90 etc)                     *
! *   - initial exclusion-list optimization (field.f90 etc)                     *
! *                                                                             *
! *******************************************************************************

!> @brief
!> - DL_MONTE-2 main program unit
!> @usage 
!> - @stdusage 
!> - @stdspecs
!> @using 
!> - `constants_module`
!> - `control_type`
!> - `control_module, only : initialise_defaults, read_control`
!> - `comms_omp_module`
!> - `comms_mpi_module`
!> - `parallel_loop_module, only : initialise_communication, initialise_loops`
!> - `field_module, only : readfield`

module DLMONTEModule

    private
    
    public :: dl_monte

    integer :: dlp_comm_world
    
    contains
    
    !> The main program of DL_MONTE-2
    subroutine dl_monte(comm_chemsh)
    
    !    use kinds_f90
        use constants_module
        use control_type
        use comms_omp_module
        use comms_mpi_module
        use control_module
        use parallel_loop_module, only : initialise_communication, initialise_loops
        use field_module, only : readfield
    
        implicit none
    
        ! ChemShell MPI communicator
    
        integer, intent(in) :: comm_chemsh
    
        type(control) :: job
    
        integer :: n_configs
        real(kind = wp) :: shortrangecut, vdwrcut, energyunit
    
        uout = 66

        ! YL: initialise ChemShell parallel communication
        call init_chemsh_comms(comm_chemsh)
        call init_omp_comms()
    
        !AB: define the only 'master' for the time being
        !AB: if replica-exchange is invoked, workgroups are created with 'master' defined per workgroup
        !master = (idnode == 0)
    
        call initialise_communication_chemshell()
    
        !open output file
        if( master ) open(uout, file='OUTPUT.000')
    
        if( master ) then
    
            write(uout,'(/,22(19x,a,/),3(19x,a,i10,a,/),36(19x,a,/))')&
                "===================================================================", &
                "* * *                                                         * * *", &
                "* * *            D L - M O N T E - v e r. 2.0 (7)             * * *", &
                "* * *                                                         * * *", &
                "* * *           STFC / CCP5 - Daresbury Laboratory            * * *", &
                "* * *                 program library  package                * * *", &
                "* * *                                                         * * *", &
                "* * *     general purpose Monte Carlo simulation software     * * *", &
                "* * *                                                         * * *", &
                "* * *     version:  2.07  / January 2020                      * * *", & !10
                "* * *     author :  J. A. Purton (2007-2015)                  * * *", &
                "* * *     website:  https://gitlab.com/dl_monte               * * *", &
                "* * *     -----------------------------------------------     * * *", &
                "* * *     co-authors:                                         * * *", &
                "* * *     -------------                                       * * *", &
                "* * *     A. V. Brukhno   - Free Energy Diff, planar pores    * * *", &
                "* * *     T. L. Underwood - Lattice/Phase-Switch MC methods   * * *", &
                "* * *     R. J. Grant     - extra potentials, optimizations   * * *", &
                "* * *     -----------------------------------------------     * * *", &
                "* * *                                                         * * *", & !20
                "* * *     current execution:                                  * * *", &
                "* * *     ==================                                  * * *", &
                "* * *     number of nodes/CPUs/cores: ",mxnode,"              * * *", &
                "* * *     number of threads per node: ",thread_num,"              * * *", &
                "* * *     number of threads employed: ",  mxthread,"              * * *", &
                "* * *                                                         * * *", &
                "* * * ======================================================= * * *", &
                "* * * Publishing data obtained with DL_MONTE? - please cite:  * * *", &
                "* * * ------------------------------------------------------- * * *", &
                "* * *                                                         * * *", & !30
                "* * *  A. V. Brukhno et al., Molecular Simulation 2019,       * * *", &
                "* * *  DOI: 10.1080/08927022.2019.1569760,                    * * *", &
                "* * *  `DL_MONTE: a multipurpose code for Monte Carlo         * * *", &
                "* * *  simulation`                                            * * *", &
                "* * *                                                         * * *", &
                "* * *  J. A. Purton, J. C. Crabtree & S. C. Parker,           * * *", &
                "* * *  Molecular Simulation 2013, 39 (14-15), 1240-1252       * * *", &
                "* * *  `DL_MONTE: A general purpose program for parallel      * * *", &
                "* * *  Monte Carlo simulation`                                * * *", &
                "* * *                                                         * * *", & !40
                "* * * ------------------------------------------------------- * * *", &
                "* * * if using Lattice-Switch method, please also cite:       * * *", &
                "* * * ------------------------------------------------------- * * *", &
                "* * *              Bruce A.D., Wilding N.B. and G.J. Ackland, * * *", &
                "* * *                         Phys. Rev. Lett. 1997, 79, 3002 * * *", &
                "* * *    `Free Energy of Crystalline Solids: A Lattice-Switch * * *", &
                "* * *                                      Monte Carlo Method`* * *", &
                "* * * ------------------------------------------------------- * * *", &
                "* * * if using Planar Pore (slit) geometry, please also cite: * * *", &
                "* * * ------------------------------------------------------- * * *", & !50
                "* * *                                         Broukhno A. V., * * *", &
                "* * *     `Free Energy and Surface Forces in polymer systems: * * *", &
                "* * *                          Monte Carlo Simulation Studies`* * *", &
                "* * *                           PhD Thesis, Lund, Sweden 2003 * * *", &
                "* * *                                                         * * *", &
                "* * *    Brukhno A. V., Akinshina A., Coldrick Z., Nelson A., * * *", &
                "* * *                                             and Auer S. * * *", &
                "* * *                              Soft Matter. 2011, 7, 1006 * * *", &
                "* * *               `Phase phenomena in supported lipid films * * *", &
                "* * *                        under varying electric potential`* * *", & !60
                "==================================================================="
                !"===================================================================", &
                !"******************************************************************"
    
        endif
    
        !get the forcefield and species
        call readfield(n_configs, shortrangecut, vdwrcut, energyunit)
    
        if (master) open(uctrl,file='CONTROL', status = 'old')
    
        call initialise_defaults(job, n_configs, shortrangecut, vdwrcut, energyunit)
        call read_control(job)
    
        if (master) close(uctrl)
    
        !setup parallel loop counters for reciprocal and two body interactions
        call initialise_loops(job%coultype, job%distribgvec, job%lparallel_atom)
    
        !start the simulation
        call run(job, n_configs)
        
        call dealloc_control(job)
    
        !AB: only master(s) close/open the files
        if( master ) then
            write(uout,*) 'normal exit '
            close(uout)
        endif
        call flush(6)
    
!        call exit_mpi_comms()
!        stop
    
    endsubroutine dl_monte
    subroutine run(job, n_confs)
    
        use constants_module,  only : uout
        use montecarlo_module, only : normalmc
        use control_type
        implicit none
    
        type(control), intent(inout) :: job
    !scp    integer, intent(in) :: n_confs : 
    !A dummy argument with the INTENT(IN) attribute shall not appear in a variable definition context 
        integer :: n_confs
    
        if( trim(adjustL(job % simulationmode)) == "trajconv" ) then
    
            !include "traj_convert.f90"
    
            call traj_convert(job, n_confs)
    
            return
    
        else if (job % useddmc) then
    
    !       call ddmc(job)
    
            call cry(uout,'(/,1x,a,/)', &
                     "ERROR: simulation mode 'DMC' is not implemented yet "//&
                    &" - nothing to do !!!",999)
    
        else
    
            call normalmc(job, n_confs)
    
        endif
    
    endsubroutine run
    subroutine init_chemsh_comms(comm_chemsh)
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    !
    ! dl_poly_3 subroutine for initialisation communication harness
    ! determining the MPI precision, and node identification and count
    !
    ! copyright - daresbury laboratory
    ! author    - i.t.todorov may 2004
    !
    !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    
!        use mpi
        use kinds_f90
        use comms_mpi_module

        Implicit None

        integer       :: comm_chemsh
        ! YL: from comms_mpi_module-PRL.f90
!        integer, save :: mpi_wp, idnode, mxnode
!        logical, save :: master, is_serial, is_parallel
!        integer       :: ierr
    
        !Call MPI_INIT(ierr)
        call MPI_COMM_DUP(comm_chemsh, dlp_comm_world, ierr)
    
        If      (wp == Selected_Real_Kind(6,30))  Then
           mpi_wp = MPI_REAL
        Else If (wp == Selected_Real_Kind(14,300))  Then
    
    ! MPI_REAL8 is apparently not in the strict MPI2 standard
    ! It is just an optional data type in the FORTRAN Bindings
           mpi_wp = MPI_DOUBLE_PRECISION
        Else If (wp == Selected_Real_Kind(30,300)) Then
           mpi_wp = MPI_REAL16
        Else
           Call error(1000)
        End If
    
        Call MPI_COMM_RANK(dlp_comm_world, idnode, ierr)
        Call MPI_COMM_SIZE(dlp_comm_world, mxnode, ierr)
    
        !AB: initially there is only one 'master'
        !AB: more are defined in parallel_loop_module::setup_workgroups()
        master = ( mxnode==1.or.idnode==0 )
    
        is_serial   = ( mxnode == 1 )
        is_parallel = .not.is_serial
    
    endsubroutine init_chemsh_comms
    !> standard/default initialisation of MPI_COMM_WORLD
    ! this routine is not called from anywhere (why ???)
    subroutine initialise_communication_chemshell()
   
!        use kinds_f90
        use parallel_loop_module
        use comms_omp_module, only : idthread, mxthread
!        use comms_mpi_module, only : master, idnode, mxnode, mpi_comm_group, mpi_comm_master
        use comms_mpi_module
   
        implicit none
   
        tb_start_atm = 1
        tb_step_atm = 1
        tb_start_mol = 1
        tb_step_mol = 1
        rcp_start = 1
        rcp_step = 1
   
        idgrp = 0
        grp_rank = 1
        wkgrp_size = 1
   
        if (mxnode > 1) then
   
            idgrp = 0
            grp_rank = idnode
            wkgrp_size = mxnode
   
            mpi_comm_group  = dlp_comm_world
            ! YL: I don't get the purpose to fake the master node using MPI_COMM_NULL. we keep it here for now 
            mpi_comm_master = MPI_COMM_NULL
   
        endif
   
        !set all processors to slave (except for single processor jobs or node 0) - this will be changed later
        master = .false.
        if (mxnode == 1 .or. idnode == 0) master = .true.
   
    endsubroutine initialise_communication_chemshell
endmodule DLMONTEModule
