polymer-nogui.c

#include <raylib.h>
#include <raymath.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>
#include <stdio.h>
#include <pthread.h>

#define SW 1500
#define SH 1000
const Vector2 mid = {.x = SW/2, .y = SH/2};
#define A(v, w) Vector2Add(v, w)
#define S(v, w) Vector2Subtract(v, w)
#define C(v, s) Vector2Scale(v, s)
#define L(v) Vector2Length(v)

#define R_C 2.5

typedef struct {
	Vector2 s, v;
} particle;
typedef unsigned int uint;

const double delta = 0.01;
double rand_uniform() { return ((double)rand() + 1.0) / ((double)RAND_MAX + 1.0); }
// Box-Muller formula
double standard_normal() { return sqrt(-2.0 * log(rand_uniform())) * cos(2.0 * M_PI * rand_uniform()); }
double   lj_force(double r, double ϵ, double σ) { return 4*ϵ/r * (12*pow(σ/r, 12) + -6*pow(σ/r, 6)); }
double ljts_force(double r, double ϵ, double σ) { return r < R_C*σ ? lj_force(r, ϵ, σ) : 0; }
double simulate_md(uint ntimes, uint N, double γ, double kₛ,
		double a, double ϵ, double σ, double kBT);
double simulate_md_vacuum(uint ntimes, uint N, double l, double E)
{ return simulate_md(ntimes, N, 1, 1000,
			l, E, l, 1000); }

int main()
{
	srand(99);
	/* InitWindow(SW, SH, "MD Simulation of Polymers"); */
	/* SetTargetFPS(60); */

	printf("eps,t,RoG,bl\n");
	for (int i = 1; i <= 10; i++) {
		simulate_md_vacuum(100000, 100, 1, i * 300);
	}

	return 0;
}

double simulate_md(uint ntimes, uint N, double γ, double kₛ,
		double a, double ϵ, double σ, double kBT)
{
	int bs = (int) ((a<1?1:a) * N);
	particle mm[N];
	for (int i = 0; i < N; i++) {
		/* mm[i].s = (Vector2) {rand()%(2*bs+1)-bs, */
		/* 	rand()%(2*bs+1)-bs}; */
		mm[i].s = (Vector2) {i*a, 0};
		mm[i].v = Vector2Zero();
	}
	double time = 0;
	/* while (!WindowShouldClose()) { */
	for (int n = 0; n < ntimes; n++) {
		for (int i = 0; i < N; i++) {
			Vector2 force = {};
			force = A(force, C(mm[i].v, -γ));
			for (int j = 0; j < N; j++) {
				if (i==j) continue;
				Vector2 r_ = S(mm[j].s, mm[i].s);
				double r = L(r_);
				if (r < 0.01) r = 0.01;
				if (abs(i-j) == 1)
					force = A(force, C(r_, 1/r * kₛ*(r-a)));
				else if (abs(i-j) >= 2)
					force = A(force, C(r_, 1/r * -ljts_force(r, ϵ, σ)));
				if (isnan(force.x))
					printf("breaking;\n");
			}
			Vector2 noise2d = {.x = standard_normal(), .y = standard_normal()};
			mm[i].v = A(mm[i].v,  A(C(force, delta),
						C(noise2d, sqrt(γ * kBT * delta))));

			// clamp down particles faster than 100
			if (L(mm[i].v) > 100) mm[i].v = C(mm[i].v, 100./L(mm[i].v));
			mm[i].s = A(mm[i].s, C(mm[i].v, delta));
		}
		Vector2 com = {};
		for (int i = 0; i < N; i++) com = A(com, C(mm[i].s, 1./N));
		for (int i = 0; i < N; i++) mm[i].s = S(mm[i].s, com);
		/* Vector2 cop = {}; */
		/* for (int i = 0; i < N; i++) cop = A(cop, C(mm[i].v, 1./N)); */
		/* for (int i = 0; i < N; i++) mm[i].v = S(mm[i].v, cop); */

		double R_g2 = 0;
		for (int i = 0; i < N; i++) R_g2 += Vector2LengthSqr(mm[i].s);
		double avgsep = 0;
		for (int i = 0; i < N-1; i++) avgsep += Vector2DistanceSqr(mm[i].s, mm[i+1].s);
		printf("%f,%f,%f,%f\n", ϵ, time, sqrt(R_g2/N), sqrt(avgsep/(N-1)));
		/* printf("%f,%f,%f,%f\n", ϵ, time, ((N-1)*a/2)/sqrt(R_g2/N), sqrt(avgsep/(N-1))/a); */

		/* BeginDrawing(); */
		/* ClearBackground(BLACK); */
		/* DrawFPS(0, 0); */
		/* for (int i = 0; i < N; i++) { */
		/* 	if (i) DrawLineV(A(mm[i-1].s, mid), A(mm[i].s, mid), GRAY); */
		/* 	DrawCircleV(A(mm[i].s, mid), 2, WHITE); */
		/* } */
		/* EndDrawing(); */

		time += delta;
	}
	/* CloseWindow(); */
	/* for (int i = 0; i < N; i++) { */
	/* 	fprintf(stderr, "[%d, %f, %f],\n", i, mm[i].s.x, mm[i].s.y); */
	/* } */
	return 0;
}