PROGRAM ideal ! demon algorithm for the one-dimensional, ideal classical gas DIM v(100) CALL initial(N,v(),mcs,E,Ed,Ecum,Edcum,accept,dvmax) FOR imcs = 1 to mcs CALL change(N,v(),E,Ed,Ecum,Edcum,accept,dvmax) NEXT imcs CALL averages(N,Ecum,Edcum,mcs,accept) END SUB initial(N,v(),mcs,E,Ed,Ecum,Edcum,accept,dvmax) RANDOMIZE INPUT prompt "number of particles = ": N INPUT prompt "initial energy of system = ": E LET Ed = 0 ! initial demon energy INPUT prompt "number of MC steps per particle = ": mcs INPUT prompt "maximum change in velocity = ": dvmax ! divide energy equally among particles LET vinitial = sqr(2*E/N) ! mass unity ! all particles have same initial velocities FOR i = 1 to N LET v(i) = vinitial NEXT i ! initialize sums LET Ecum = 0 LET Edcum = 0 LET accept = 0 END SUB SUB change(N,v(),E,Ed,Ecum,Edcum,accept,dvmax) FOR i = 1 to N LET dv = (2*rnd - 1)*dvmax ! trial change in velocity LET ipart = int(rnd*N + 1) ! select random particle LET vtrial = v(ipart) + dv ! trial velocity ! trial energy change LET de = 0.5*(vtrial*vtrial - v(ipart)*v(ipart)) IF de <= Ed then LET v(ipart) = vtrial LET accept = accept + 1 LET Ed = Ed - de LET E = E + de END IF NEXT i ! accumulate data after each Monte Carlo step per particle LET Ecum = Ecum + E LET Edcum = Edcum + Ed END SUB SUB averages(N,Ecum,Edcum,mcs,accept) LET norm = 1/mcs LET Edave = Edcum*norm ! mean demon energy LET norm = norm/N LET accept_prob = accept*norm ! acceptance probability ! system averages per particle LET Esave = Ecum*norm ! mean energy per system particle PRINT "mean demon energy ="; Edave PRINT "mean system energy per particle ="; Esave PRINT "acceptance probability ="; accept_prob END SUB