sim1 = function(n=40, mu=0.6343) {
  a  = rnorm(n,   0.0, 1)
  b1 = rnorm(n/2, mu - mu/4, 1)
  b2 = rnorm(n/2, mu + mu/4, 1) # interaction is mu/2

  dat = data.frame(Group = c(rep("C", n), rep("T", n)),
                   Sex   = rep(c(rep("F", n/2), rep("M", n/2)), 2),
                   Val   = c(a, b1, b2))
  fit   = lm(Val ~ Group*Sex, dat)
  return(summary(fit)$coefficients)
}

sim2 = function(n=40, mu=0.6343) {
  a  = rnorm(n,   0.0, 1)
  b1 = rnorm(n/2, mu - mu/4, 1)
  b2 = rnorm(n/2, mu + mu/4, 1) # interaction is mu/2

  dat = data.frame(Group = c(rep("C", n), rep("T", n)),
                   Sex   = rep(c(rep("M", n/2), rep("F", n/2)), 2),  # we only exchange labels M & F
                   Val   = c(a, b1, b2))
  fit   = lm(Val ~ Group*Sex, dat)
  return(summary(fit)$coefficients)
}

n = 40
d = power.t.test(n=n, power=0.8)$delta

# variation 1
res1=sapply(1:1000, function(x) sim1(n, d))
mean(res1[14,]<0.05) # power main effect

# variation 2
res2=sapply(1:1000, function(x) sim2(n, d))
mean(res2[14,]<0.05) # power main effect