I ran a protein Fep simulation of Valine to alanine bound to a protein using 30 windows. When I analyzed the simulation, I got high peak of lambda windows between 0.96 and 1
Attached below is the plot when I added 1 and without 1

Also below is my code

lam_vals = [round(i / 29, 4) for i in range(30)] # [0.0000, 0.0625, ..., 1.0000]

Loop over all lambda values

for lam in lam_vals:
lam_dir = f"lambda_{lam:.4f}" # e.g., lambda_0.0000
print(f"\n🔁 Running lambda = {lam:.4f}")
# --- Minimisation ---
min_protocol = BSS.Protocol.FreeEnergyMinimisation(lam=lam, lam_vals=lam_vals, steps=1000)
min_process = BSS.Process.Gromacs(system=solvated_system_holo, protocol=min_protocol, work_dir=f"{lam_dir}/min")
min_process.start()
min_process.wait()
minimised_system = min_process.getSystem()
# --- Equilibration ---
eq_protocol = BSS.Protocol.FreeEnergyEquilibration(lam=lam, lam_vals=lam_vals, runtime="125ps")
eq_process = BSS.Process.Gromacs(system=minimised_system, protocol=eq_protocol, work_dir=f"{lam_dir}/eq")
eq_process.start()
eq_process.wait()
equilibrated_system = eq_process.getSystem()
# --- Production ---
prod_protocol = BSS.Protocol.FreeEnergyProduction(lam=lam, lam_vals=lam_vals, runtime="4ns")
prod_process = BSS.Process.Gromacs(system=equilibrated_system, protocol=prod_protocol, work_dir=f"{lam_dir}/prod")
prod_process.start()
prod_process.wait()
print(f"✅ Finished lambda = {lam:.4f}")