anchor.diff:                  m%xx%ind(w) = xx_cur
                   m%yy%ind(w) = yy_cur
                   m%zz%ind(w) = zz_cur
  -                call random_move(m,w,xx_cur,yy_cur,zz_cur,xx_new,yy_new,zz_new, max_move)
  +                call random_move(m, w, xx_cur, yy_cur, zz_cur, xx_new, yy_new, zz_new, max_move)
               end do
               ! Update hutches, data for chi2, and chi2/del_chi
  -            call hutch_move_atom(m,w,xx_new, yy_new, zz_new)
  +            call hutch_move_atom(m, w, xx_new, yy_new, zz_new)
       
               call eam_mc(m, w, xx_cur, yy_cur, zz_cur, xx_new, yy_new, zz_new, te2)
               te2 = te2/m%natoms
  -            !if(myid .eq. 0) write(*,*) "Energy = ", te2
   
               ! Calculate a randnum for accept/reject
               randnum = ran2(iseed2)
  -            ! Decide whether to reject just based on the energy
  -            accepted = .true.
  -            if(chi2_initial*0.1 > chi2_no_energy) then
  -            if(log(1.0-randnum) > -(te2-chi2_old)*beta) then
  -                accepted = .false.
  -                call reject_position(m, w, xx_cur, yy_cur, zz_cur)
  -                call hutch_move_atom(m,w,xx_cur, yy_cur, zz_cur)  !update hutches.  
  -                e2 = e1 ! eam
  -                if(myid.eq.0) write(*,*) "MC move rejected solely due to energy."
  -            endif
  -            endif
  -                
  -            if(accepted) then
  +
               call fem_update(m, w, res, k, vk, v_background, scatfact_e, mpi_comm_world, istat, square_pixel)
  -            !write(*,*) "I am core", myid, "and I have exited from fem_update into the main hrmc block."
   
               chi2_no_energy = chi_square(alpha, vk_exp, vk_exp_err, vk, scale_fac, nk)
  -
  -
               chi2_new = chi2_no_energy + te2
  -            !chi2_new = chi2_new * m%natoms
               del_chi = chi2_new - chi2_old
               call mpi_bcast(del_chi, 1, mpi_double, 0, mpi_comm_world, mpierr)
   
  
positive.code:         ! Write initial vk 
          open(unit=52,file=trim(vki_fn),form='formatted',status='unknown')
              do i=1, nk
                  write(52,*) k(i),vk(i)
              enddo
          close(52)
      endif
  
      !------------------- Call potential energy function. -----------------!
  
      if(use_hrmc) then ! End here if we only want femsim. Set the variable above.
  
      call read_eam(m,eam_filename)
      call eam_initial(m,te1)
      te1 = te1/m%natoms
  
negative.code:             write(*,*) "   Temperature =      ", temperature
              write(*,*) "   Max Move=          ", max_move
              write(*,*)
              ! Reset energy/chi_squared file
              open(36,file=trim(chi_squared_file), form='formatted', status='unknown')
                  write(36,*) "step, chi2, energy"
                  write(36,*) i, chi2_no_energy, te1
              close(36)
              open(37,file=trim(acceptance_rate_fn), form='formatted', status='unknown', access='append')
                  write(37,*) "step, acceptance rate averaged over last 1000 steps"
              close(37)
          endif
  