anchor.diff:  #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
  +    write(*,*) "Successfully initialized MPI"
  +    communicator = mpi_comm_world
  +    call mpi_comm_rank(communicator, myid, mpierr)
  +    call mpi_comm_size(communicator, numprocs, mpierr)
  +    write(*,*) "I am rank", myid, "of", numprocs
  +
  +    call mpi_barrier(communicator, mpierr)
   #else
       myid = 0
  -    colored_comm = 0
  +    communicator = 0
       numprocs = 1
   #endif
   
  
positive.code:             endif
  
              if(energy_accepted) then
                  call fem_update(m, w, res, k, vk, scatfact_e, colored_comm, istat)
  
                  chi2_no_energy = chi_square(alpha, vk_exp, vk_exp_err, vk, scale_fac, nk)
                  chi2_new = chi2_no_energy + te2
                  del_chi = chi2_new - chi2_old
  #ifndef SERIAL
                  call mpi_bcast(del_chi, 1, mpi_double, 0, colored_comm, mpierr)
  #endif
  
                  ! Test if the move should be accepted or rejected based on del_chi
  
negative.code:     !------------------- Read inputs and initialize. -----------------!
  
      ! Set input filenames.
      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)
  