path: root/common.h

/* See copyright information at the end of the file */

/*! @file common.h
 * @brief Common definitions needed in simulation code.
 */

#ifndef __COMMON_H__
#define __COMMON_H__

enum mode {
	MODE_UNKNOWN = 0,
	MODE_ANGLES = 1 << 1,
	MODE_TRANSLATE = 1 << 2,
};

enum log_level {
	NO_LOG,
	WARN,
	LOG,
	VERBOSE,
};

static enum log_level log_level = WARN;
#ifdef DEBUG
#define DBG(fmt, ...)                                                        \
	do {                                                                 \
		if (log_level > NO_LOG)                                      \
			fprintf(stderr, "%s:%d:0 [WARN] " fmt, __FILE__,     \
				__LINE__, __VA_ARGS__);                      \
	} while (0)
#define LOG(fmt, ...)                                                        \
	do {                                                                 \
		if (log_level > WARN)                                        \
			fprintf(stderr, "%s:%d:0 [LOG] " fmt, __FILE__,      \
				__LINE__, __VA_ARGS__);                      \
	} while (0)
#define VERBOSE(fmt, ...)                                                    \
	do {                                                                 \
		if (log_level > LOG) {                                       \
			fprintf(stderr, "%s:%d:0 [VERBOSE] " fmt, __FILE__,  \
				__LINE__, __VA_ARGS__);                      \
		}                                                            \
	} while (0)
#else
#define DBG(fmt, ...)                                                        \
	do {                                                                 \
		if (log_level > NO_LOG)                                      \
			fprintf(stderr, "[WARNING] " fmt, __VA_ARGS__);      \
	} while (0)
#define LOG(fmt, ...)                                                        \
	do {                                                                 \
		if (log_level > WARN)                                        \
			fprintf(stderr, "[LOG] " fmt, __VA_ARGS__);          \
	} while (0)
#define VERBOSE(fmt, ...)                                                    \
	do {                                                                 \
		if (log_level > LOG) {                                       \
			fprintf(stderr, "[VERBOSE] " fmt, __VA_ARGS__);      \
		}                                                            \
	} while (0)
#endif /* DEBUG */

/*
 * NOTE: nsdirname = directory name without trailing slash
 *       sdirname = directory name with trailing slash
 */

#define MIN(a, b) ((a) < (b) ? (a) : (b))
#define EPS 1e-10
#define STRLEN(x) (sizeof(x) / sizeof(char) - 1)
#define STEP 0.1
#define METADATA_FILENAME "METADATA"
#define METADATA_FILENAME_LEN STRLEN(METADATA_FILENAME)
#define ANGLES_FILENAME "angles.energies"
#define ANGLES_FILENAME_LEN STRLEN(ANGLES_FILENAME)
#define NORMALIZED_FILENAME "normalized.energies"
#define NORMALIZED_FILENAME_LEN STRLEN(NORMALIZED_FILENAME)
#define NORMALIZED_MOVABLE_FILENAME "normalized_movable.energies"
#define NORMALIZED_MOVABLE_FILENAME_LEN STRLEN(NORMALIZED_MOVABLE_FILENAME)
#define NON_NORMALIZED_FILENAME "non-normalized.energies"
#define NON_NORMALIZED_FILENAME_LEN STRLEN(NON_NORMALIZED_FILENAME)
#define RIGID_FILENAME "rigid.energies"
#define RIGID_FILENAME_LEN STRLEN(RIGID_FILENAME)
#define LATTICE_FILENAME "lattice_rXXX.in"
#define LATTICE_FILENAME_LEN STRLEN(LATTICE_FILENAME)
#define RUN_FILENAME "run.sh"
#define RUN_FILENAME_LEN STRLEN(RUN_FILENAME)
#define RIGID_INPUT_FILENAME "rigid.in"
#define RIGID_INPUT_FILENAME_LEN STRLEN(RIGID_INPUT_FILENAME)
#define TRANSLATION_FILENAME "translation.energies"
#define TRANSLATION_FILENAME_LEN STRLEN(TRANSLATION_FILENAME)
#define GEOMETRIC_BETA 0.85
#define EXPONENTIAL_BETA 0.85
#define LUNDY_MEES_BETA 0.005
#define LINEAR_NAME "linear"
#define GEOMETRIC_NAME "geometric"
#define LUNDY_MEES_NAME "lundy-mees"
#define EXPONENTIAL_NAME "exponential"

typedef void (*CoolingFn)(FILE *out, int ntry, int damping, int seed,
	int runtime, double temperature, int annealing, double beta);

typedef struct Metadata Metadata;
struct Metadata {
	int annealing; /*!< Number of consecutive simulated annealing steps */
	int damping;   /*!< Default langevin damping
			  (https://doc.lammps.org/fix_langevin.html) */
	int from;      /*!< First twist angle in rotation simulation */
	int natoms;    /*!< Number of movable atoms in the sample */
	int nstep;     /*!< Number of angular steps for "angular" simulation,
			* number of translational steps in "translate"
			* simulation.
			*/
	int ntry;    /*!< Number of times the simulation should be repeated */
	int runtime; /*!< Cooling time for simulated annealing */
	int seed;    /*!< Default langevin seed */
	int to;	     /*!< Final twist angle in rotation simulation */
	int tot_natoms;	 /*!< Total number of atoms in the sample */
	size_t timestep; /*!< Simulation timestep in LAMMPS units */
	double angle; /*!< Initial rotation angle in "translate" simulation */
	double apot;  /*!< Characteristic length scale of the quasiperiodic
			 potential */
	double acol;  /*!< Characteristic length scale of colloidal bonds (in
			 um) */
	double beta;  /*!< Cooling schedule parameter, see each implementation
		       * for the  default value
		       */
	double bond_coeff; /*!< K Harmonic coefficient used in the bond_coeff
			    * command
			    * (https://doc.lammps.org/bond_harmonic.html) */
	double step; /*!< Angular step in rotation simulations, i.e. the total
			number of angular configurations is (to - from)/step
		      */
	double temperature;   /*!< Maximum temperature for simulated annealing
			       */
	double tolerance;     /*!< Threshold value to stop minimization */
	double v0;	      /*!< Default potential in fJ */
	CoolingFn cooling_fn; /*!< One of linear, geometric, lundy mees or
				 exponential  */
	String min_style;     /*!< One of fire, cg, sd, quickmin */
	String data;	      /*!< Path to LAMMPS initial data file */
	double translation[2]; /*!< Center of the simulation in units of apot
				*/
	enum mode mode; /*!< Bitmask of the possible simulation modes */
};

/*! Fill out *mode, *v0, *natoms and *nstep read from
 * nsdirname/METADATA_FILENAME.
 *
 * @param[in] r Global region.
 * @param[in] nsdirname Simulation directory (without trailing slash).
 * @param[out] m Result of metadata parsing.
 * @return true if no errors occurred , false otherwise.
 */
bool dir_metadata(Region *r, const char nsdirname[static restrict 1],
	Metadata m[static restrict 1]);
/*! Parse metadata string.
 * Expected metadata format is space-separated pair on each line e.g.:
 * @code
 * MODE RIGID
 * ATOMS 9407
 * V0 1e-8
 * NSTEP 60
 * ...
 * @endcode
 *
 * @param[in] r Memory region for allocations.
 * @param[in] metadata Contents of metadata file.
 * @param[in] nsdirname Simulation directory (without trailing slash).
 * @param[in] path Name of the metadata file.
 * @param[out] m Result of the parsing.
 * @return true if no errors occurred , false otherwise.
 */
bool get_metadata(Region *r, String metadata,
	const char nsdirname[static restrict 1],
	const char path[static restrict 1], Metadata m[static restrict 1]);
/*! Strip '/' from the end of sdirname. */
void strip_slash(char sdirname[static 1]);
/*! Read the LAMMPS data file filename, filling *natoms with the
 * number of movable atoms and *tot_natoms with the total number of
 * atoms, returning true if no errors are encountered.
 */
bool extract_natoms(Region *r, const char filename[static 1],
	int natoms[static 1], int tot_natoms[static 1]);
/*! Write metadata information.
 *
 * The format is the following:
 *
 * @code
 * MODE ANGLES
 * V0 1e-8
 * ANNEALING 1
 * ATOMS 9407
 * ...
 * @endcode
 *
 * @param[in] r Global Region.
 * @param[in] nsdirname Simulation directory (without trailing slash).
 * @param[in] data LAMMPS data file.
 * @param[in] m Metadata that needs to be recorded in nsdirname.
 * @return true if no errors occurred, false otherwise.
 */
bool write_metadata(Region *r, const char nsdirname[static 1],
	const char data[static 1], Metadata m[static 1]);
/*! Return the name corresponding to the specified cooling schedule.
 * @param[in] fn Cooling function.
 * @return Cooling function name.
 */
const char *cooling_schedule_name(CoolingFn fn);
/*! Return pointer to cooling function from its name.
 * @param[in] name The cooling function's name.
 * @return Pointer to the cooling function.
 */
CoolingFn cooling_schedule_from_name(const char *name);
/*! \f$T_k = exp(-\frac{\beta*k}{tot})*T_0\f$ */
void cooling_exponential(FILE out[static 1], int ntry, int damping, int seed,
	int runtime, double temperature, int annealing, double beta);
/*! \f$T_k = \beta^k*T_0\f$ */
void cooling_geometric(FILE out[static 1], int ntry, int damping, int seed,
	int runtime, double temperature, int annealing, double beta);
/*! LAMMPS default when running the langevin thermostat */
void cooling_linear(FILE out[static 1], int ntry, int damping, int seed,
	int runtime, double temperature, int annealing, double beta);
/*! \f$T_{k+1} = \frac{T_k}{1+\beta*T_k}\f$ */
void cooling_lundy_mees(FILE out[static 1], int ntry, int damping, int seed,
	int runtime, double temperature, int annealing, double beta);

#ifdef COMMON_IMPLEMENTATION
void
strip_slash(char sdirname[static 1])
{
	size_t dirlen = strlen(sdirname);
	if (sdirname[dirlen - 1] == '/') {
		sdirname[dirlen - 1] = '\0';
	}
}

bool
dir_metadata(Region *r, const char nsdirname[static restrict 1],
	Metadata m[static restrict 1])
{
	FILE *fp;
	char *metadata;
	size_t metadata_len;
	String contents = { 0 };

	metadata_len = strlen(nsdirname) + METADATA_FILENAME_LEN + 2;
	metadata = region_malloc(r, metadata_len);
	snprintf(metadata, metadata_len, "%s/" METADATA_FILENAME, nsdirname);

	if (!(fp = fopen(metadata, "r"))) {
		DBG("Error opening %s: %s\n", metadata, strerror(errno));
		goto err;
	}

	if (!string_read(r, fp, &contents)) {
		DBG("Error reading from %s: %s\n", metadata, strerror(errno));
		goto err;
	}

	LOG("Reading metadata from %s\n", metadata);
	if (!get_metadata(r, contents, nsdirname, metadata, m)) {
		DBG("Error getting metadata from: %s\n", metadata);
		goto err;
	}

	fclose(fp);
	return true;
err:
	if (fp)
		fclose(fp);
	return false;
}

bool
get_metadata(Region *r, String metadata,
	const char nsdirname[static restrict 1],
	const char path[static restrict 1], Metadata m[static restrict 1])
{
	String line = { 0 };

	m->mode = MODE_UNKNOWN;
	m->nstep = m->natoms = m->tot_natoms = -1;
	m->step = STEP;

	line = string_next_line(&metadata);
	while (line.str) {
		String tok;
		tok = string_next_tok(&line, ' ');
		string_trim(&tok);
		if (string_cmp(tok, STRING_LITERAL("MODE")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected mode, found nothing.\n",
					path);
				return false;
			}
			string_trim(&tok);
			if (string_cmp(tok, STRING_LITERAL("ANGLES")) == 0) {
				m->mode = MODE_ANGLES;
				VERBOSE("Mode: " String_Fmt "\n",
					String_Arg(tok));
			} else if (string_cmp(
					   tok, STRING_LITERAL("TRANSLATE"))
				== 0) {
				m->mode = MODE_TRANSLATE;
				VERBOSE("Mode: " String_Fmt "\n",
					String_Arg(tok));
			} else {
				DBG("%s: Unknown mode: " String_Fmt
				    ", skipping\n",
					path, String_Arg(tok));
				return false;
			}
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("ATOMS")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected number of atoms, found "
				    "nothing.\n",
					path);
				return false;
			}
			VERBOSE("atoms: " String_Fmt "\n", String_Arg(tok));
			string_trim(&tok);
			if (!string_strtoi(tok, &m->natoms)) {
				DBG("%s: Expected number of atoms, "
				    "found: " String_Fmt "\n",
					path, String_Arg(tok));
				return false;
			}
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("V0")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected potential, found "
				    "nothing.\n",
					path);
				return false;
			}
			VERBOSE("v0: " String_Fmt "\n", String_Arg(tok));
			string_trim(&tok);
			if (!string_strtod(tok, &m->v0)) {
				DBG("%s: Expected v0, "
				    "found: " String_Fmt "\n",
					path, String_Arg(tok));
				return false;
			}
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("TEMPERATURE")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected temperature, found "
				    "nothing.\n",
					path);
				goto next;
			}
			VERBOSE("temperature: " String_Fmt "\n",
				String_Arg(tok));
			string_trim(&tok);
			if (!string_strtod(tok, &m->temperature)) {
				DBG("%s: Expected temperature, "
				    "found: " String_Fmt "\n",
					path, String_Arg(tok));
			}
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("APOT")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected apot, found "
				    "nothing.\n",
					path);
				goto next;
			}
			VERBOSE("apot: " String_Fmt "\n", String_Arg(tok));
			string_trim(&tok);
			if (!string_strtod(tok, &m->apot)) {
				DBG("%s: Expected apot, "
				    "found: " String_Fmt "\n",
					path, String_Arg(tok));
			}
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("ACOL")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected acol, found "
				    "nothing.\n",
					path);
				goto next;
			}
			VERBOSE("acol: " String_Fmt "\n", String_Arg(tok));
			string_trim(&tok);
			if (!string_strtod(tok, &m->acol)) {
				DBG("%s: Expected acol, "
				    "found: " String_Fmt "\n",
					path, String_Arg(tok));
			}
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("BOND_COEFFICIENT"))
			== 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected bond coefficient, found "
				    "nothing.\n",
					path);
				goto next;
			}
			VERBOSE("bond_coeff: " String_Fmt "\n",
				String_Arg(tok));
			string_trim(&tok);
			if (!string_strtod(tok, &m->bond_coeff)) {
				DBG("%s: Expected bond_coefficient, "
				    "found: " String_Fmt "\n",
					path, String_Arg(tok));
			}
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("MIN-STYLE")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected minimization style, found "
				    "nothing.\n",
					path);
				goto next;
			}
			VERBOSE("min_style: " String_Fmt "\n",
				String_Arg(tok));
			string_trim(&tok);
			m->min_style = tok;
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("NSTEP")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected number of steps, found "
				    "nothing.\n",
					path);
				return false;
			}
			VERBOSE("nstep: " String_Fmt "\n", String_Arg(tok));
			string_trim(&tok);
			if (!string_strtoi(tok, &m->nstep)) {
				DBG("%s: Expected nstep, found: " String_Fmt
				    "\n",
					path, String_Arg(tok));
				return false;
			}
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("TOT_NATOMS")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected total number of atoms, "
				    "found "
				    "nothing.\n",
					path);
				goto next;
			}
			VERBOSE("tot_natoms: " String_Fmt "\n",
				String_Arg(tok));
			string_trim(&tok);
			if (!string_strtoi(tok, &m->tot_natoms)) {
				DBG("%s: Expected total number of atoms, "
				    "found: " String_Fmt "\n",
					path, String_Arg(tok));
			}
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("ANNEALING")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected number of annealing steps, "
				    "found"
				    "nothing.\n",
					path);
				goto next;
			}
			VERBOSE("annealing: " String_Fmt "\n",
				String_Arg(tok));
			string_trim(&tok);
			if (!string_strtoi(tok, &m->annealing)) {
				DBG("%s: Expected number of annealing steps, "
				    "found: " String_Fmt "\n",
					path, String_Arg(tok));
			}
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("STEP")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected angular step size, found "
				    "nothing.\n",
					path);
				return false;
			}
			VERBOSE("step: " String_Fmt "\n", String_Arg(tok));
			string_trim(&tok);
			if (!string_strtod(tok, &m->step)) {
				DBG("%s: Expected angular step size, "
				    "found: " String_Fmt "\n",
					path, String_Arg(tok));
				return false;
			}
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("COOLING_FN")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected cooling function, found "
				    "nothing.\n",
					path);
				goto next;
			}
			VERBOSE("cooling_fn: " String_Fmt "\n",
				String_Arg(tok));
			string_trim(&tok);
			m->cooling_fn = cooling_schedule_from_name(tok.str);
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("TX")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected x component of translation "
				    "vector, found nothing.\n",
					path);
				goto next;
			}
			VERBOSE("tx: " String_Fmt "\n", String_Arg(tok));
			string_trim(&tok);
			if (!string_strtod(tok, &m->translation[0])) {
				DBG("%s: Expected x component of translation "
				    "vector, found: " String_Fmt "\n",
					path, String_Arg(tok));
				return false;
			}
			goto next;
		}
		if (string_cmp(tok, STRING_LITERAL("TY")) == 0) {
			tok = string_next_tok(&line, ' ');
			if (!tok.str) {
				DBG("%s: Expected y component of translation "
				    "vector, found nothing.\n",
					path);
				goto next;
			}
			VERBOSE("ty: " String_Fmt "\n", String_Arg(tok));
			string_trim(&tok);
			if (!string_strtod(tok, &m->translation[1])) {
				DBG("%s: Expected y component of translation "
				    "vector, found: " String_Fmt "\n",
					path, String_Arg(tok));
				return false;
			}
			goto next;
		}

		String tok2 = string_next_tok(&line, ' ');
		if (!tok2.str) {
			DBG("%s: Non terminated metadata entry: " String_Fmt
			    ".\n",
				path, String_Arg(tok));
			return false;
		}

	next:
		line = string_next_line(&metadata);
	}
	if (m->tot_natoms == -1 || m->tot_natoms == m->natoms) {
		/* This means we used the old version of simulation.c
		 * to generate this simulation, we need to update the
		 * metadata accordingly: so we must overwrite the old
		 * METADATA.
		 */
		DIR *d = opendir(nsdirname);
		struct dirent *dir = NULL;
		size_t nslen = strlen(nsdirname);
		char *data = NULL, *metadata_filename = NULL;
		int errnum = 0;
		if (!d) {
			return NULL;
		}
		VERBOSE("Number of movable atoms not found%s\n", "");
		while ((dir = readdir(d))) {
			if (dir->d_type == DT_REG) {
				String sname = make_string(
					strlen(dir->d_name), dir->d_name);
				if (string_ends_with(
					    STRING_LITERAL(".lmpdat"),
					    sname)) {
					data = region_malloc(r,
						(nslen + 1 + sname.count + 1)
							* sizeof(char));
					memset(data, 0,
						nslen + 1 + sname.count + 1);
					strncpy(data, nsdirname, nslen);
					data[nslen] = '/';
					strncpy(&data[nslen + 1], dir->d_name,
						sname.count);
					break;
				}
			}
		}
		VERBOSE("Parsing %s to find the number of fixed atoms.\n",
			data);
		if (!write_metadata(r, "/tmp", data, m)) {
			DBG("Error updating %s/" METADATA_FILENAME ": %s.\n",
				nsdirname, strerror(errno));
			return false;
		}
		VERBOSE("tot_natoms: %d\n", m->tot_natoms);
		VERBOSE("natoms: %d\n", m->natoms);
		metadata_filename = region_malloc(
			r, nslen + 1 + METADATA_FILENAME_LEN + 1);
		memset(metadata_filename, 0,
			nslen + 1 + METADATA_FILENAME_LEN + 1);
		memcpy(metadata_filename, nsdirname, nslen);
		metadata_filename[nslen] = '/';
		memcpy(&metadata_filename[nslen + 1], METADATA_FILENAME,
			METADATA_FILENAME_LEN);
		/* TODO(mario): This is the last thing I was working
		 * on... I'm tired and it's late.
		 *
		 * I still need to look at a couple of things:
		 *
		 * 1. Why does the data filename appear in the
		 * METADATA file?
		 *
		 * 2. Why does analyze complain about not being able
		 * to parse the metadata? I should look at
		 * get_metadata or dir_metadata for this
		 */
		if ((errnum = cp(
			     "/tmp/" METADATA_FILENAME, metadata_filename))
			!= 0) {
			DBG("Error copying /tmp/" METADATA_FILENAME
			    " to %s: %s. Continuing anyway.\n",
				metadata_filename, strerror(errnum));
		}
		/* Recursive calls suck here, but it is what it is. I
		 * need to find a way to modify tot_natoms here...
		 * Maybe I should extract_natoms here and pass natoms
		 * and tot_natoms to write_metadata, but then I would
		 * need to do it everywhere */
		/* return get_metadata(r, metadata, nsdirname, path, m); */
	}

	return m->mode != MODE_UNKNOWN && m->natoms != -1 && m->nstep != -1
		&& m->tot_natoms != -1;
}

bool
extract_natoms(Region *r, const char filename[static 1], int natoms[static 1],
	int tot_natoms[static 1])
{
	String contents = { 0 }, line = { 0 };
	bool in_atoms = false;
	int fixed_atoms = 0;

	if (!string_read_file(r, filename, &contents)) {
		DBG("Error reading file %s: %s\n", filename, strerror(errno));
		return false;
	}

	line = string_next_line(&contents);
	while (line.str) {
		if (string_ends_with(STRING_LITERAL(" atoms"), line)) {
			String tok = { 0 };

			tok = string_next_tok(&line, ' ');
			if (!string_strtoi(tok, natoms)) {
				DBG("Expected number of atoms in %s, "
				    "found: " String_Fmt "\n",
					filename, String_Arg(tok));
				return false;
			}
		} else if (string_starts_with(
				   STRING_LITERAL("Atoms"), line)) {
			in_atoms = true;
			string_next_line(&contents); /* Skip empty line. */
		} else if (in_atoms) {
			String tok = { 0 };
			int atom_type = 0;
			if (!line.str
				|| string_starts_with(
					STRING_LITERAL("Bonds"), line)) {
				*tot_natoms = *natoms;
				*natoms -= fixed_atoms;
				return true;
			} else {
				tok = string_next_tok(&line, ' ');
				tok = string_next_tok(&line, ' ');
				tok = string_next_tok(&line, ' ');
				if (!string_strtoi(tok, &atom_type)) {
					DBG("Expected a number as the bond "
					    "type, found " String_Fmt "\n",
						String_Arg(tok));
					return false;
				}
				if (atom_type == 2) {
					fixed_atoms++;
				}
			}
		}

		line = string_next_line(&contents);
	}

	DBG("Error: could not find number of atoms in %s.\n", filename);
	return false;
}

bool
write_metadata(Region *r, const char nsdirname[static 1],
	const char data[static 1], Metadata m[static restrict 1])
{
	char *metadata_path;
	size_t metadata_len, dirlen;
	FILE *metadata = NULL;

	dirlen = strlen(nsdirname);
	metadata_len = dirlen + METADATA_FILENAME_LEN + 2;
	metadata_path = region_malloc(r, metadata_len);
	if (snprintf(metadata_path, metadata_len, "%s/" METADATA_FILENAME,
		    nsdirname)
		< 0) {
		DBG("Encoding error.%s", "\n");
		goto err;
	}

	if (!(metadata = fopen(metadata_path, "w"))) {
		DBG("Error opening %s file for writing: %s\n", metadata_path,
			strerror(errno));
		goto err;
	}

	if (!extract_natoms(r, data, &m->natoms, &m->tot_natoms)) {
		DBG("Error calculating number of atoms from data %s\n", data);
		goto err;
	}

	fprintf(metadata, "MODE ");
	switch (m->mode) {
	case MODE_ANGLES:
		fprintf(metadata, "ANGLES\n");
		break;
	case MODE_TRANSLATE:
		fprintf(metadata, "TRANSLATE\n");
		break;
	default:
		assert(0 && "Unreachable");
	}
	fprintf(metadata,
		"V0 %.6e\n"
		"ATOMS %d\n"
		"TOT_NATOMS %d\n"
		"TEMPERATURE %f\n"
		"ANNEALING %d\n"
		"NSTEP %d\n"
		"MIN-STYLE " String_Fmt "\n" /* Need to upcase */
		"BOND-COEFFICIENT %.6e\n"
		"STEP %.2e\n"
		"APOT %.4e\n"
		"ACOLL %.4e\n"
		"TX %.4e\n"
		"TY %.4e\n",
		m->v0, m->natoms, m->tot_natoms, m->temperature, m->annealing,
		m->nstep, String_Arg(m->min_style), m->bond_coeff, m->step,
		m->apot, m->acol, m->translation[0], m->translation[1]);
	if (m->cooling_fn == cooling_exponential) {
		fprintf(metadata, "BETA . %.3f\n",
			m->beta > 0 ? m->beta : EXPONENTIAL_BETA);
	} else if (m->cooling_fn == cooling_lundy_mees) {
		fprintf(metadata, "BETA . %.3f\n",
			m->beta > 0 ? m->beta : LUNDY_MEES_BETA);
	} else if (m->cooling_fn == cooling_geometric) {
		fprintf(metadata, "BETA . %.3f\n",
			m->beta > 0 ? m->beta : GEOMETRIC_BETA);
	}

	fclose(metadata);
	return true;
err:
	if (metadata)
		fclose(metadata);
	return false;
}

const char *
cooling_schedule_name(CoolingFn fn)
{
	if (fn == cooling_linear) {
		return LINEAR_NAME;
	} else if (fn == cooling_geometric) {
		return GEOMETRIC_NAME;
	} else if (fn == cooling_exponential) {
		return EXPONENTIAL_NAME;
	} else if (fn == cooling_lundy_mees) {
		return LUNDY_MEES_NAME;
	} else {
		return "unknown";
	}
}

CoolingFn
cooling_schedule_from_name(const char *name)
{
	if (strcmp(name, LINEAR_NAME) == 0) {
		return cooling_linear;
	} else if (strcmp(name, GEOMETRIC_NAME) == 0) {
		return cooling_geometric;
	} else if (strcmp(name, LUNDY_MEES_NAME) == 0) {
		return cooling_lundy_mees;
	} else if (strcmp(name, EXPONENTIAL_NAME) == 0) {
		return cooling_exponential;
	}
	return NULL;
}

void
cooling_exponential(FILE out[static 1], int ntry, int damping, int seed,
	int runtime, double temperature, int annealing, double beta)
{
	int random = rand();

	if (ntry > 1) {
		fprintf(out,
			"variable initial_temperature equal %.6e/10/${try}\n",
			temperature);
	} else {
		fprintf(out, "variable initial_temperature equal %.6e\n",
			temperature);
	}
	fprintf(out,
		"variable beta equal %.3f\n"
		"variable last equal %d\n"
		"variable nannealing equal %d\n"
		"variable loop_%d loop 1 ${nannealing}\n"
		"label label_%d\n",
		beta > 0 ? beta : EXPONENTIAL_BETA, annealing - 1, annealing,
		random, random);

	/* Start of Loop */
	fprintf(out,
		"variable high_temp equal "
		"exp(-${beta}*(${loop_%d}-1)/"
		"${nannealing})*${initial_temperature}\n"
		"if \"${loop_%d} < ${last}\" then &\n"
		"    \"variable low_temp equal "
		"exp(-${beta}*${loop_%d}/"
		"${nannealing})*${initial_temperature}\" &\n"
		"else &\n"
		"    \"variable low_temp equal 0\"\n"
		"fix annealing interno langevin ${high_temp} ${low_temp} %d "
		"%d\n"
		"fix e interno nve\n"
		"fix 2d interno enforce2d\n"
		"run %d\n"
		"unfix annealing\n"
		"unfix e\n"
		"unfix 2d\n"
		"variable high_temp delete\n"
		"variable low_temp delete \n",
		random, random, random, damping, seed, runtime / annealing);

	/* End of loop */
	fprintf(out,
		"next loop_%d\n"
		"jump SELF label_%d\n",
		random, random);
}

void
cooling_geometric(FILE out[static 1], int ntry, int damping, int seed,
	int runtime, double temperature, int annealing, double beta)
{
	int random = rand();

	if (ntry > 1) {
		fprintf(out,
			"variable initial_temperature equal %.6e/10/${try}\n",
			temperature);
	} else {
		fprintf(out, "variable initial_temperature equal %.6e\n",
			temperature);
	}
	fprintf(out,
		"variable beta equal %.3f\n"
		"variable last equal %d\n"
		"variable loop_%d loop 1 %d\n"
		"label label_%d\n",
		beta > 0 ? beta : GEOMETRIC_BETA, annealing - 1, random,
		annealing, random);

	/* Start of Loop */
	fprintf(out,
		"variable high_temp equal "
		"${beta}^${loop_%d}*${initial_temperature}\n"
		"if \"${loop_%d} < ${last}\" then &\n"
		"    \"variable low_temp equal "
		"${beta}^(1+${loop_%d})*${initial_temperature}\" &\n"
		"else &\n"
		"    \"variable low_temp equal 0\"\n"
		"fix annealing interno langevin ${high_temp} ${low_temp} %d "
		"%d\n"
		"fix e interno nve\n"
		"fix 2d interno enforce2d\n"
		"run %d\n"
		"unfix annealing\n"
		"unfix e\n"
		"unfix 2d\n"
		"variable high_temp delete\n"
		"variable low_temp delete\n",
		random, random, random, damping, seed, runtime / annealing);

	/* End of loop */
	fprintf(out,
		"next loop_%d\n"
		"jump SELF label_%d\n",
		random, random);
}

void
cooling_linear(FILE out[static 1], int ntry, int damping, int seed,
	int runtime, double temperature, int annealing, double beta)
{
	if (ntry > 1) {
		fprintf(out,
			"variable initial_temperature equal %.6e/10/${try}\n",
			temperature);
	} else {
		fprintf(out, "variable initial_temperature equal %.6e\n",
			temperature);
	}
	fprintf(out,
		"fix annealing interno langevin ${initial_temperature} 0 %d "
		"%d\n"
		"fix e interno nve\n"
		"fix 2d interno enforce2d\n"
		"run %d\n",
		damping, seed, runtime / annealing);
}

void
cooling_lundy_mees(FILE out[static 1], int ntry, int damping, int seed,
	int runtime, double temperature, int annealing, double beta)
{
	int random = rand();

	if (ntry > 1) {
		fprintf(out,
			"variable initial_temperature equal %.6e/10/${try}\n",
			temperature);
	} else {
		fprintf(out, "variable initial_temperature equal %.6e\n",
			temperature);
	}
	fprintf(out,
		"variable beta equal %.3f\n"
		"variable last equal %d\n"
		"variable loop_%d loop 1 %d\n"
		"label label_%d\n",
		beta > 0 ? beta : LUNDY_MEES_BETA, annealing - 1, random,
		annealing, random);

	/* Start of Loop */
	fprintf(out,
		"if \"${loop_%d} < 2\" then &\n"
		"\"variable low_temp equal ${initial_temperature}\" &\n"
		"    \"variable high_temp equal ${initial_temperature}\" "
		"&\n"
		"else &\n"
		"    \"variable high_temp equal ${low_temp}\"\n"
		"variable low_temp delete\n"
		"variable low_temp equal "
		"${high_temp}/(1+${beta}*${high_temp})\n"
		"fix annealing interno langevin ${high_temp} ${low_temp} %d "
		"%d\n"
		"fix e interno nve\n"
		"fix 2d interno enforce2d\n"
		"run %d\n"
		"unfix annealing\n"
		"unfix e\n"
		"unfix 2d\n"
		"variable high_temp delete\n",
		random, damping, seed, runtime / annealing);

	/* End of loop */
	fprintf(out,
		"next loop_%d\n"
		"jump SELF label_%d\n",
		random, random);
}

#endif /* COMMON_IMPLEMENTATION */
#endif /* __COMMON_H__ */

/*
 * Copyright ©️ 2023 Mario Forzanini <mf@marioforzanini.com>
 *
 * This file is part of my bachelor thesis.
 *
 * This file is free software: you can redistribute it and/or modify it
 * under the terms of the GNU General Public License as published by the
 * Free Software Foundation, either version 3 of the License, or (at your
 * option) any later version.
 *
 * This file is distributed in the hope that it will be useful, but
 * WITHOUT ANY WARRANTY; without even the implied warranty of
 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
 * General Public License for more details.
 *
 * You should have received a copy of the GNU General Public License
 * along with this file. If not, see <https://www.gnu.org/licenses/>.
 *
 */