anchor.diff:      !------------------- Read inputs and initialize. -----------------!
   
       ! Set input filenames.
  -#ifndef SERIAL
  -    call mpi_bcast(param_filename, 256, MPI_CHARACTER, 0, colored_comm, mpierr)
  -#endif
       if(myid.eq.0) write(*,*) "Paramfile: ", trim(param_filename)
  +#ifndef SERIAL
  +    call mpi_bcast(param_filename, 256, MPI_CHARACTER, 0, communicator, mpierr)
  +#endif
       
       ! Read input parameters
       call read_inputs(param_filename, femsim, model_filename, femfile, Q, ntheta, nphi, npsi, scale_fac_initial, eam_filename, step_start, step_end, temp_move_decrement, temperature, max_move, cutoff_r, seed, alpha, vk_exp, k, vk_exp_err)
  
positive.code: #endif
  
  #ifndef SERIAL
      ! Read in this special environemnt variable set by OpenMPI
      ! Context:  http://stackoverflow.com/questions/35924226/openmpi-mpmd-get-communication-size
      ! (See the solution by Hristo Iliev)
      ! The below two lines set `color` to be 0, 1, ..., P where P is the number
      ! of executables given on the command line within `mpiexec` that are colon-deliminated
      call get_environment_variable("OMPI_MCA_orte_app_num", color_str)
      call str2int(color_str, color, istat)
  
      call mpi_init_thread(MPI_THREAD_MULTIPLE, ipvd, mpierr)
      call mpi_comm_rank(mpi_comm_world, myid, mpierr)
      call mpi_comm_size(mpi_comm_world, numprocs, mpierr)
  
      ! Split the world communicator into separate pieces for each mpiexec subprogram
      call mpi_barrier(mpi_comm_world, mpierr)
      call mpi_comm_split(mpi_comm_world, color, myid, colored_comm, mpierr)
  
      ! Get each core's rank within its new communicator
      call mpi_comm_rank(colored_comm, myid, mpierr)
      ! Get the total number of cores (ie the sum of all the worlds)
      call mpi_comm_size(colored_comm, numprocs, mpierr)
  
      call mpi_barrier(colored_comm, mpierr)
      write(*,'(A10, I2, A4, I2, A11, I2, A6, I2, A18, I2)') "I am core ", myid, " of ", numprocs, " with color", color, ", root", 0, ", and communicator", colored_comm
  #else
      myid = 0
      colored_comm = 0
      numprocs = 1
  #endif
  
  
negative.code: 
  #ifndef SERIAL
      ! Free the sub-communicators and finalize
      call mpi_comm_get_parent(communicator, mpierr)
      if(communicator .ne. mpi_comm_null) then
          call mpi_comm_disconnect(communicator, mpierr)
          write(*,*) "Disconnected from parent!"
      endif
      call mpi_finalize(mpierr)
  #endif
  
  end program hrmc
  
  