Initial commit.

Simple simulation using the verlet integration with electrostatic non-bound potentials

Author: TerminalCursor

LICENSE

@@ -0,0 +1,26 @@
+Copyright (c) <year>, <copyright holder>
+All rights reserved.
+
+Redistribution and use in source and binary forms, with or without
+modification, are permitted provided that the following conditions are met:
+
+1. Redistributions of source code must retain the above copyright notice, this
+   list of conditions and the following disclaimer.
+2. Redistributions in binary form must reproduce the above copyright notice,
+   this list of conditions and the following disclaimer in the documentation
+   and/or other materials provided with the distribution.
+
+THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS IS" AND
+ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED
+WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE
+DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER OR CONTRIBUTORS BE LIABLE FOR
+ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES
+(INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
+LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND
+ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
+(INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
+SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
+
+The views and conclusions contained in the software and documentation are those
+of the authors and should not be interpreted as representing official policies,
+either expressed or implied, of the FreeBSD Project.

Makefile

@@ -0,0 +1,25 @@
+# $@ = target file
+# $< = first dependency
+# $^ = all dependencies
+#
+C_SOURCES = $(wildcard Physics/*.c *.c)
+C_HEADERS = $(wildcard Physics/*.h *.h)
+C_OBJECTS = ${C_SOURCES:.c=.o}
+
+CC = gcc
+CFLAGS=-lm
+
+# FIRST RULE IS RUN BY DEFAULT
+build: ${C_OBJECTS}
+	${CC} -o $@ $^ ${CFLAGS}
+
+%.o: %.c ${C_HEADERS}
+	${CC} -c $< -o $@
+
+clean:
+	rm -rf ${C_OBJECTS} build
+
+backup:
+	cp -r ./ $(shell pwd)_bak
+
+.PHONY: backup clean

Physics/Differentiation.c

@@ -0,0 +1,48 @@
+#include "Differentiation.h"
+
+Particle** copyPpp(int n, Particle** p) {
+	Particle** copy = (Particle**)malloc(sizeof(Particle*) * n);
+	for(int i = 0; i < n; i++) {
+		copy[i] = malloc(sizeof(Particle));
+		copy[i] = copy_particle(p[i]);
+	}
+	return copy;
+}
+
+void clonePpp(int n, Particle** src, Particle** dest) {
+	for(int i = 0; i < n; i++) {
+		for(int j = 0; j < DIM; j++) {
+			dest[i]->pos[j] = src[i]->pos[j];
+			dest[i]->vel[j] = src[i]->vel[j];
+		}
+		dest[i]->m = src[i]->m;
+		dest[i]->q = src[i]->q;
+	}
+}
+
+void freePpp(int n, Particle** p) {
+	for(int i = 0; i < n; i++) {
+		free_particle(p[i]);
+	}
+	free(p);
+}
+
+double* particle_force(int n, Particle** fin, int j, double (*potential)(int, Particle**, int)) {
+	double* returnSlope = calloc(DIM, sizeof(double));
+	for(int i = 0; i < n; i++) {
+		if(j != i) {
+			for(int k = 0; k < DIM; k++) {
+				Particle** pCopyUp = copyPpp(n, fin);
+				Particle** pCopyDown = copyPpp(n, fin);
+				pCopyUp[i]->pos[k] += DERIV_STEP;
+				pCopyDown[i]->pos[k] -= DERIV_STEP;
+				double high = potential(n, pCopyUp, j);
+				double low = potential(n, pCopyDown, j);
+				freePpp(n, pCopyUp);
+				freePpp(n, pCopyDown);
+				returnSlope[k] += (high - low) / (2 * DERIV_STEP);
+			}
+		}
+	}
+	return returnSlope;
+}

Physics/Differentiation.h

@@ -0,0 +1,18 @@
+#ifndef DIFFERENTIATION_H
+#define DIFFERENTIATION_H
+
+#include <stdlib.h>
+#include <math.h>
+#include "Particle.h"
+
+#define DERIV_STEP 0.0000001
+
+double* particle_force(int, Particle**, int, double (*potential)(int, Particle**, int));
+
+Particle** copyPpp(int n, Particle** p);
+
+void clonePpp(int n, Particle** src, Particle** dest);
+
+void freePpp(int n, Particle** p);
+
+#endif

Physics/Particle.c

@@ -0,0 +1,29 @@
+#include "Particle.h"
+
+Particle* make_particle_at(int n, double* pos) {
+	Particle* p = malloc(sizeof(Particle));
+	for(int i = 0; i < n; i++) {
+		if(i < DIM) {
+			p->pos[i] = pos[i];
+			p->vel[i] = 0;
+		}
+	}
+	p->m = 1;
+	p->q = 1;
+	return p;
+}
+
+Particle* copy_particle(Particle* p) {
+	Particle* copy = malloc(sizeof(Particle));
+	for(int i = 0; i < DIM; i++) {
+		copy->pos[i] = p->pos[i];
+		copy->vel[i] = p->vel[i];
+	}
+	copy->m = p->m;
+	copy->q = p->q;
+	return copy;
+}
+
+void free_particle(Particle* p) {
+	free(p);
+}

Physics/Particle.h

@@ -0,0 +1,20 @@
+#ifndef PARTICLE_H
+#define PARTICLE_H
+
+#include <stdlib.h>
+
+#define DIM 3
+
+typedef struct {
+	double pos[DIM];
+	double vel[DIM];
+	double m, q;
+} Particle;
+
+Particle* make_particle_at(int, double*);
+
+Particle* copy_particle(Particle*);
+
+void free_particle(Particle*);
+
+#endif

Physics/Potentials.c

@@ -0,0 +1,16 @@
+#include "Potentials.h"
+
+double electrostaticPotential(int n, Particle** ps, int j) {
+	double retPot = 0;
+	for(int i = 0; i < n; i++) {
+		if(i != j) {
+			double rs = 0;
+			for(int k = 0; k < DIM; k++) {
+				rs += pow(ps[i]->pos[k] - ps[j]->pos[k], 2);
+			}
+			rs = sqrt(rs);
+			retPot += ps[i]->q * ps[j]->q * 9e9 / rs;
+		}
+	}
+	return retPot;
+}

Physics/Potentials.h

@@ -0,0 +1,9 @@
+#ifndef POTENTIALS_C
+#define POTENTIALS_C
+
+#include <math.h>
+#include "Particle.h"
+
+double electrostaticPotential(int n, Particle** ps, int j);
+
+#endif

Physics/Verlet.c

@@ -0,0 +1,44 @@
+#include "Verlet.h"
+
+Particle* next_position(int n, Particle** prev, Particle** current, int i, int nf, double (*potential[])(int, Particle**, int)) {
+	Particle* next = copy_particle(current[i]);
+	// From the verlet algorithm - first terms for the next position
+	for(int k = 0; k < DIM; k++) {
+		next->pos[k] = 2 * current[i]->pos[k] - prev[i]->pos[k];
+	}
+	// Calculate the forces
+	for(int j = 0; j < n; j++) {
+		if(i != j) {
+			for(int k = 0; k < nf; k++) {
+				double* force = particle_force(n, current, i, potential[k]);
+				for(int m = 0; m < DIM; m++) {
+					next->pos[m] += force[m] / next->m * DT * DT;
+				}
+				free(force);
+			}
+		}
+	}
+	return next;
+}
+
+Particle* next_velocity(int n, Particle** prev, Particle** next, int i) {
+	Particle* updated = copy_particle(next[i]);
+	for(int k = 0; k < DIM; k++) {
+		updated->vel[k] = (1 / (2.0 * DT)) * (updated->pos[k] - prev[i]->pos[k]);
+	}
+	return updated;
+}
+
+void next_step(int n, Particle** prev, Particle** current, Particle** next, int nf, double (*potentials[])(int, Particle**, int)) {
+	Particle** next_p = (Particle**)malloc(sizeof(Particle*) * n);
+	for(int i = 0; i < n; i++) {
+		next_p[i] = next_position(n, prev, current, i, nf, potentials);
+	}
+	for(int i = 0; i < n; i++) {
+		next[i] = next_velocity(n, prev, next_p, i);
+	}
+	for(int i = 0; i < n; i++) {
+		free_particle(next_p[i]);
+	}
+	free(next_p);
+}

Physics/Verlet.h

@@ -0,0 +1,17 @@
+#ifndef VERLET_H
+#define VERLET_H
+
+#include <math.h>
+#include <stdlib.h>
+#include "Particle.h"
+#include "Differentiation.h"
+
+#define DT 0.00001
+
+Particle* next_position(int n, Particle** prev, Particle** current, int i, int nf, double (*potential[])(int, Particle**, int));
+
+Particle* next_velocity(int n, Particle** prev, Particle** next, int i);
+
+void next_step(int n, Particle** prev, Particle** current, Particle** next, int nf, double (*potential[])(int, Particle**, int));
+
+#endif

main.c

@@ -0,0 +1,56 @@
+#include <stdio.h>
+#include <math.h>
+#include "Physics/Differentiation.h"
+#include "Physics/Particle.h"
+#include "Physics/Verlet.h"
+#include "Physics/Potentials.h"
+
+double mean_dist(int n, Particle** ps) {
+	double mean = 0;
+	for(int i = 0; i < n; i++) {
+		for(int j = i+1; j < n; j++) {
+			double rs = 0;
+			for(int k = 0; k < DIM; k++) {
+				rs += pow(ps[i]->pos[k]-ps[j]->pos[k], 2);
+			}
+			rs = sqrt(rs);
+			mean += rs;
+		}
+	}
+	return mean;
+}
+
+int main(int argc, const char* argv[]) {
+	printf("Booting up a simple simulation that only considers the non-bonded electrostatic potential\n");
+	int MAX_STEP = 10;
+	double (*pf[1])(int, Particle**, int);
+	pf[0] = electrostaticPotential;
+	double pos[] = { 0, 0, 0 };
+	int len = 2;
+	Particle* p1 = make_particle_at(3, pos);
+	p1->q = -1;
+	pos[0] = 2;
+	Particle* p2 = make_particle_at(3, pos);
+	Particle* rs[] = { p1, p2 };
+
+	Particle** prev = copyPpp(len, rs);
+	Particle** current = copyPpp(len, rs);
+	Particle** next = copyPpp(len, rs);
+	for(int step = 0; step < MAX_STEP; step++) {
+		for(int i = 0; i < len; i++) {
+			printf("%i's position is now: (%f, %f, %f)\n", i, next[i]->pos[0], next[i]->pos[1], next[i]->pos[2]);
+			printf("%i's velocity is now: (%f, %f, %f)\n", i, next[i]->vel[0], next[i]->vel[1], next[i]->vel[2]);
+		}
+		printf("Mean distance is: %f\n", mean_dist(len, next));
+		clonePpp(len, current, prev);
+		clonePpp(len, next, current);
+		next_step(len, prev, current, next, 1, pf);
+	}
+	freePpp(len, next);
+	freePpp(len, current);
+	freePpp(len, prev);
+	for(int i = 0; i < len; i++) {
+		free_particle(rs[i]);
+	}
+	return 0;
+}