solve.py

# translated from https://gitlab.com/manzhara157/nanobelts/-/blob/main/example.pas

import array
import math
import random
import os
from contextlib import suppress


n0 = 200000

z = 12

tmax = 1000000
dt = 0.1

ra_kinet = 1
rb_kinet = 1
rc_kinet = 1

ra_therm = 1
rb_therm = 1
rc_therm = 1


ra_bal = 1
rb_bal = 1
rc_bal = 1

j = 0.008
xeq = 0.01
d0 = 0.01

a0 = 200.0
b0 = 200.0
c0 = 200.0
vtot = n0 * a0 * b0 * c0 * 100

lich0 = 100

t = 0
d = d0

ra_bal = 1 / 2
rb_bal = math.sqrt(2)
rc_bal = math.sqrt(2)

sumv = 0
sump = 0

state = [0] * (n0 + 1)
a = [0] * (n0 + 1)
b = [0] * (n0 + 1)
c = [0] * (n0 + 1)

with open("a.txt", "r") as af:
    with open("b.txt", "r") as bf:
        with open("c.txt", "r") as cf:
            for i in range(1, n0 + 1):
                state[i] = 1
                a[i] = float(af.readline())
                b[i] = float(bf.readline())
                c[i] = float(cf.readline())
                sumv = sumv + a[i] * b[i] * c[i]
                sump = sump + ((a[i] * b[i]) + (b[i] * c[i]) + (c[i] * a[i]))


v0 = sumv
n = n0

lich = lich0


def clean_files():
    with suppress(FileNotFoundError):
        os.remove("mean_abc_py.txt")
    with suppress(FileNotFoundError):
        os.remove("fmean_v_py.txt")
    with suppress(FileNotFoundError):
        os.remove("fd_py.txt")
    with suppress(FileNotFoundError):
        os.remove("fn_py.txt")
    with suppress(FileNotFoundError):
        os.remove("fp_py.txt")
    with suppress(FileNotFoundError):
        os.remove("fv_py.txt")


def save_data():
    ln_t = math.log(t)
    with open("mean_abc_py.txt", "a") as f:
        f.write(
            "%f %f %f %f\n"
            % (ln_t, math.log(mean_a), math.log(mean_b), math.log(mean_c))
        )

    with open("fmean_v_py.txt", "a") as f:
        f.write("%f %f\n" % (ln_t, math.log(mean_v)))

    with open("fd_py.txt", "a") as f:
        f.write("%f %f\n" % (ln_t, math.log(d)))

    with open("fn_py.txt", "a") as f:
        f.write("%f %f\n" % (ln_t, math.log(n)))

    with open("fp_py.txt", "a") as f:
        f.write("%f %f\n" % (ln_t, math.log(p)))

    with open("fv_py.txt", "a") as f:
        f.write("%f %f\n" % (ln_t, math.log(sumv)))


clean_files()

random.seed(31337)

while t < tmax:
    t = t + dt
    suma = 0
    sumb = 0
    sumc = 0
    sumv = 0
    sump = 0

    neg = 0

    i = 1

    anew = [0] * (n + 1)
    bnew = [0] * (n + 1)
    cnew = [0] * (n + 1)

    while i <= n:
        anew[i] = a[i] + dt * ra_kinet * (
            d - (rc_therm / c[i]) - (rb_therm / b[i]) - (ra_bal / ra_kinet) * j
        )
        bnew[i] = b[i] + dt * rb_kinet * (
            d - (ra_therm / a[i]) - (rc_therm / c[i]) - (rb_bal / rb_kinet) * j
        )
        cnew[i] = c[i] + dt * rc_kinet * (
            d - (ra_therm / a[i]) - (rb_therm / b[i]) - (rc_bal / rc_kinet) * j
        )

        if (anew[i] < 0) or (bnew[i] < 0) or (cnew[i] < 0):
            a[i] = a[n]
            b[i] = b[n]
            c[i] = c[n]
            n = n - 1
            neg += 1
        else:
            suma = suma + anew[i]
            sumb = sumb + bnew[i]
            sumc = sumc + cnew[i]
            sumv = sumv + anew[i] * bnew[i] * cnew[i]
            sump = sump + (
                (anew[i] * bnew[i]) + (bnew[i] * cnew[i]) + (cnew[i] * anew[i])
            )
            i = i + 1

    a = anew
    b = bnew
    c = cnew

    d = ((xeq + d0) * (1 - v0 / vtot) + v0 / vtot - sumv / vtot) / (
        1 - sumv / vtot
    ) - xeq
    p = 2 * sump

    if lich == lich0:
        mean_a = suma / n
        mean_b = sumb / n
        mean_c = sumc / n
        mean_v = sumv / n

        save_data()

        lich = 0

    lich = lich + 1