# LAMMPS Script for Dumbbell with Random State Transitions

# Variables
variable T equal 300
variable total_steps equal 50000    # Total simulation steps
variable current_step equal 0       # Counter for the current simulation step

# Transition probabilities
variable p_AtoB equal 0.01  # Probability to transition from A to B

variable p_BtoA equal 0.01  # Probability to transition from B to A



# Set initial state
variable current_state equal 0  # Start in state A (0)
variable next_state equal 0

# Initialization
units lj
dimension 3
boundary p p p
atom_style bond
timestep 0.000001

# Read the data file
read_data dumbell_input.dat

# Define masses
mass 1 0.000001
mass 2 0.000001

# Create groups for different atom types
group type1 type 1
group type2 type 2

# Define bond style and initial coefficients
bond_style harmonic
bond_coeff 1 500.0 0.5  # Initial bond length for state A

# Langevin dynamics for type1 and type2
fix langevin1 type1 brownian 0.0001 12908410 gamma_t 1 rng gaussian
fix langevin2 type2 brownian 0.0001 12908410 gamma_t 10 rng gaussian  # Initial damping for state A

# Compute the bond distance
compute bond_dist all bond/local dist

# Dump bond distances and positions
dump 1 all local 1 bond_dist.dat c_bond_dist
dump 2 all custom 100 dump.dumbbell id type x y z
dump 3 all xyz 10 dumbbellxyz.xyz

# Create a group for both particles
group all_particles union type1 type2

# Compute the center of mass
compute com_all all_particles com
thermo 100
thermo_style custom step c_com_all[1] c_com_all[2] c_com_all[3]

# Simulation loop
label loop

# Generate a random number
variable rand equal random(0,1,12345)

# Transition logic for state A (0)
if "${current_state} == 0" then "variable next_state equal 0"
if "${current_state} == 0 && ${rand} < ${p_AtoB}" then "variable next_state equal 1"


# Transition logic for state B (1)
if "${current_state} == 1" then "variable next_state equal 1"
if "${current_state} == 1 && ${rand} < ${p_BtoA}" then "variable next_state equal 0"



# Update the current state
variable current_state equal ${next_state}

# Print current state
print "Current step: ${current_step}, Current state: ${current_state}"

# Update bond length and damping based on the new state
if "${current_state} == 0" then "bond_coeff 1 500.0 0.5"
if "${current_state} == 1" then "bond_coeff 1 500.0 2.0"


unfix langevin2
if "${current_state} == 0" then "fix langevin2 type2 brownian 0.0001 12908410 gamma_t 10 rng gaussian"
if "${current_state} == 1" then "fix langevin2 type2 brownian 0.0001 12908410 gamma_t 50 rng gaussian"


# Run one timestep
run 1

# Increment step counter
variable current_step equal ${current_step}+1

# Check if total simulation time is reached
if "${current_step} >= ${total_steps}" then "jump SELF end"

# Loop back
jump SELF loop

label end
print "Simulation completed after ${total_steps} steps."