////////////////////////////////////////////////////////////////////////
#ifndef NETWORK_HPP
#define NETWORK_HPP

#include <iostream>
#include <iomanip>
#include <fstream>
#include <vector>
#include <thread>
#include <string>
#include <sstream>
#include <mutex>

#define _USE_MATH_DEFINES
#include <math.h>

#include <stdio.h>
// #include <intrin.h>

#include "helper.hpp"
#include "randomnumbers.h"

using namespace std;

////////////////////////////////////////////////////////////////////////
std::mutex mtx;
// for random number generation
long idum;				// for the seed of the random numbers
float ran2(long *idum); // generates random number in [0,1)

// https://helloacm.com/the-rdtsc-performance-timer-written-in-c/
//  Windows
#ifdef _WIN32
#include <intrin.h>
// ''Read Time-Stamp Counter'' returns the number of clock cycles since last reset
unsigned long long rdtsc()
{
	return __rdtsc();
}
//  Linux/GCC
#else
unsigned long long rdtsc()
{
	unsigned int lo, hi;
	__asm__ __volatile__("rdtsc"
						 : "=a"(lo), "=d"(hi));

	return ((unsigned long)hi << 32) | lo;
}

#endif

/* #include <x86intrin.h>
// ''Read Time-Stamp Counter'' returns the number of clock cycles since last reset
unsigned long long rdtsc()
{
	return __rdtscp();
} */

/* #ifdef _WIN32
#include <intrin.h>
#else
#include <x86intrin.h>
#endif
// ''Read Time-Stamp Counter'' returns the number of clock cycles since last reset
unsigned long long rdtsc()
{
	return __rdtsc();
} */

////////////////////////////////////////////////////////////////////////
class network
{
private:
	const double lc = 7.0;								  // crystal length
	const size_t init_dist = 7;							  // initial distance between nearest neighbor nodes /  "lattice constant"
	const double d_init = static_cast<double>(init_dist); // double value of lattice constant
	const double nodes_rcut = 0.5 * d_init * 2.41;		  // "cut-off radius" for connecting nodes
	const double nodes_rcutsq = nodes_rcut * nodes_rcut;  // cut-off is estimated by first and second nearest neighbor distances:
														  // fist nearest neighbor distance is a = d_init
														  // second nearest neighbor distance is b = d_init * sqrt(2), with sqrt(2) \approx 1.41
														  // -> cut-off estimation is 0.5 * (a + b)

	const double r0 = 0.5 * d_init; // characteristic distance for interaction
	const double r0_inv = 1.0 / r0;

	const double theta0 = 1.3; // free energy for crystallization
	const double theta0_c = theta0 * lc;

	size_t NodesN;		   // initial number of nodes (approximated),
	size_t Links_Node;	   // number of links per node estimated by dimension
	double Links_per_Node; // actual number of links per node after setting up the initial lattice
	size_t dim;			   // dimension of the network

	vector<double> box_size;					// vector containing the lengths of the box
	double init_length;							// initial length of the square/cubic simulation box
	vector<vector<double>> Nodes;				// vector containing the positions of the Nodes
	vector<vector<double>> Segments;			// vector containing the number of segments in a link
	vector<vector<size_t>> Neighbors;			// vector containing the indices of nearest-neighbor nodes
	vector<vector<vector<double>>> Periodic_BC; // vector containing the information whether a link crosses the boundary of the simulation box
	vector<vector<bool>> Defects;				// vector containing the information whether a link is crystallizable
	bool crystallizing;							// information whether the network is able to crystallize
	vector<double> Crystals;					// vector containing the number of crystals in a link

	vector<vector<size_t>> Connections; // vector containing the indices of the two nodes which are connected by a link
										// and the type of the link, i.e. 0=polymer-polymer, 1=filler-filler, 2=polymer-filler

	double kappa;		// measure on how inhomogenous the initial lattice is
	double defect_frac; // fraction of non-crystallizable links

	bool filled;
	double filler_frac;	 // fraction of nodes considered as filler
	size_t filler_bonds; // number of filler bonds

	vector<vector<size_t>> Link_Neighbors; // vector containing the indices of the neighboring links for each link
	vector<vector<double>> Link_Centers;   // vector containing the positions of the center of ech link
	size_t LinksN;						   // total number of links
	double rcut;						   // cut-off radius for the interaction of the links
	double rcutsq;
	double rm; // "skin" radius for neighbor list of links
	double rmsq;

	double max_strain;	 // maximum strain
	double strain_steps; // steps in which the strain is adjusted

	bool make_force_hist;
	vector<double> force_hist_lambda; // strains at which the forces on each link are saved

	bool local_link_data_save; // stretch, crystallinity, extension and energy of each link will be saved; NOT for pf and ff bonds
	vector<vector<double>> local_link_data;

	bool Save_Segments; // decide to save number of segments in each link

	vector<vector<double>> Force; // vector containing the force acting on each node
	vector<vector<double>> Force_abs;
	double G;					// total free energy
	vector<double> Free_Energy; // vector containing the free energy for each node
	double eta;					// enthalpic interaction
	double eta_c;

	double A0_inv;	   // inverse of original cross section of the sample
	double ratio_A1A2; // ratio of the cross sections A1 and A2 dependent on lambda
	double sigma;	   // shear stress in stretching direction

	double chi; // crystallinity index

	bool RuptureCheck;		// allow links to crack
	double RuptureEnergy;	// critical energy per Kuhn segment for rupture of a link
	double RuptureCryst;	// maximum crystallinity that is allowed for links which break at the critical force
	bool CrackCentered;		// whether the crack is centered in the network or at the boundary
	bool CrackCheck;		// whether a crack is inserted initially in the network
	double CrackSize;		// crack size in units of approximately twice the initial distance of the nodes
	size_t LinksCracked;	// number of links that broke in a stretching step
	vector<bool> CrackArea; // vector containing the information whether a node is in the area around the inserted crack

	vector<bool> LinkCheck; // an element of this vector is true if the length of a link r is larger than the number n of its segments
	size_t LinksCrackedRN;	// number of links with r>n that cracked
	fstream LinksRuptured;	// save n, r, c of the links that ruptured

	// parameters of the FIRE algorithm for energy minimization
	const double FIRE_alpha_start = 0.95;
	const double FIRE_dt_start = 0.01;
	const double FIRE_dt_max = 1.0;
	const double FIRE_f_inc = 1.1;
	const double FIRE_f_dec = 0.9;
	const double FIRE_f_alpha = 0.9;
	const double FIRE_mass = 10.0;
	const size_t FIRE_N_min = 10;
	const size_t FIRE_N_max = 10000;
	const double FIRE_w_rel = pow(10.0, -8.0);
	const double FIRE_mass_inv = 1.0 / FIRE_mass;

	// for saving of the data to plot the network
	bool save;
	// size_t save_iteration;
	size_t save_counter;
	fstream savePos;
	fstream saveCryst;
	fstream saveBox;

	bool read_initial_config;

	// for morphology generator
	bool use_morphology_generator;
	vector<bool> Filler;							// vector containing information whether a node is considered as filler or not
	size_t MCattempts;								// number of MC interchange attempts per node for the calculation of the total number of MC steps
	size_t MCsteps;									// number of MC steps for morphology generation
	const double surface_tension_polymer_d = 20.24; // for NR in mJ/m^2
	const double surface_tension_polymer_p = 5.46;
	const double surface_tension_filler_d = 27.0; // 18.7;	//27.0
	const double surface_tension_filler_p = 0.0;  // 22.7;	//0.0
	const double surface_tension_polymer = surface_tension_polymer_d + surface_tension_polymer_p;
	const double surface_tension_filler = surface_tension_filler_d + surface_tension_filler_p;
	const double interface_tension_polymer_filler = surface_tension_polymer + surface_tension_filler - 2.0 * (std::sqrt(surface_tension_filler_d * surface_tension_polymer_d) + std::sqrt(surface_tension_filler_p * surface_tension_polymer_p));
	const double interface_area = 0.421; // = a/(k_B*T) in m^2/mJ
	const double interface_polymer_filler_factor = interface_tension_polymer_filler * interface_area;
	vector<vector<size_t>> filler_neighbors;
	const double r_morphology = 2.0 * d_init; //
	const double r_morphology_sq = r_morphology * r_morphology;

	const double spring_ff = 5.0;	   // spring constant for filler-filler interaction
	const double spring_pf = 4.0;	   // spring constant for polymer-filler interaction
	const double spring_pf_weak = 1.0; // spring constant for polymer-filler interaction at large range
	const bool cut_ff = true;
	const bool cut_pf = true;
	const double spring_ff_half = 0.5 * spring_ff;
	const double spring_pf_half = 0.5 * spring_pf;
	const double spring_pf_weak_half = 0.5 * spring_pf_weak;
	const double spring_ff_inv = 3.0 / spring_ff; // for computation of equilibrium distances
	const double spring_pf_inv = 3.0 / spring_pf;
	size_t no_ff; // total number of filler-filler bonds
	size_t no_pf; // total number of polymer-filler bonds

	const double R_factor_ff = 1.1; // factors which are multiplied with the equilibrium distances to obtain cut in the potential
	const double R_factor_pf = 1.1;
	const double R_factor_ff_sq = R_factor_ff * R_factor_ff;
	const double R_factor_pf_sq = R_factor_pf * R_factor_pf;
	vector<vector<double>> eq_dist; // vector containing the equilibrium distances of the links, R*r_0 and whether the nodes are currently connected or not
									// (reversible bond breaking for filler-filler and polymer-filler links)

	const double filler_rcut = nodes_rcut; // r_morphology; // cut-off radius for interaction of filler nodes
	const double filler_rcutsq = filler_rcut * filler_rcut;

	double sigma1;
	double sigma2;

	////////////////////////////////////////////////////////////////////////

	// this function reads the initial positions of the nodes from a file
	// The data in this file has to be of the same form as in the files in which the positons are saved
	// if the morphology generator is used. The positions have to correspond to them after the first energy minimization.
	// The linkages are set up in "read_initial_links()", where also filler nodes are defined.
	void read_initial_node_positions(const std::string &_initial_positions_file)
	{
		std::cout << "# Reading " + _initial_positions_file + "..." << std::endl;

		std::ifstream file(_initial_positions_file);
		std::string line;
		std::vector<double> matrix;

		if (!file)
		{
			std::cout << "# Cannot find " + _initial_positions_file + "." << std::endl;
			exit(0);
		}

		while (std::getline(file, line))
		{
			std::stringstream ss(line);
			double line_data;

			double value;
			while (ss >> value)
			{
				line_data = value;
			}
			matrix.push_back(line_data);
		}
		file.close();

		if (matrix.size() != dim * NodesN)
		{
			std::cout << "# Number of positions given in the file does not fit the given parameters." << std::endl;
			exit(0);
		}

		// set the positions of the nodes
		for (size_t i = 0; i < dim; i++)
		{
			for (size_t j = 0; j < NodesN; j++)
			{
				const size_t index = (i * NodesN) + j;
				Nodes[j].push_back(matrix[index]);
			}
		}
		matrix.clear();
		std::cout << "# Reading " + _initial_positions_file + " successful." << std::endl;
	}

	// this function reads the initial links from a file and assigns the filler type to certain nodes
	// The equilibrium distances of filler-filler and polymer-filler bonds are also determined.
	// They are given by the current distances of the nodes.
	void read_initial_links(const std::string &_initial_linkage_file)
	{
		std::cout << "# Reading " + _initial_linkage_file + " ..." << std::endl;

		std::ifstream file(_initial_linkage_file);
		std::string line;
		std::vector<std::vector<double>> matrix;

		if (!file)
		{
			std::cout << "# Cannot find " + _initial_linkage_file << "." << std::endl;
			std::cin.ignore();
			exit(0);
		}

		while (std::getline(file, line))
		{
			std::stringstream ss(line);
			std::vector<double> line_data;

			double value;
			while (ss >> value)
			{
				line_data.push_back(value);
			}
			matrix.push_back(line_data);
		}
		file.close();

		// set up the links
		Neighbors.resize(NodesN);
		Segments.resize(NodesN);
		Filler.resize(NodesN);
		Defects.resize(NodesN);
		LinksN = matrix.size();
		Crystals.clear();
		Connections.clear();
		Periodic_BC.resize(NodesN);
		for (size_t idx = 0; idx < LinksN; idx++)
		{
			const size_t j = matrix[idx][0];
			const size_t k = matrix[idx][1];
			const size_t n = matrix[idx][2];
			const size_t defect_link = matrix[idx][3];
			bool defect_info = false;
			if (crystallizing == true && defect_link == 1)
			{
				defect_info = true;
			}

			const bool filler_j = matrix[idx][4];
			const bool filler_k = matrix[idx][5];

			Filler[j]=filler_j;
			Filler[k]=filler_k;
			
			// save node k as a neighbor of node j
			Neighbors[j].push_back(k);
			// save number of Kuhn segments in the link connecting node j and k
			Segments[j].push_back(n);
			// save information whether the considered link is a "defect"
			Defects[j].push_back(defect_info);
			// set number of crystals in the considered link to zero
			Crystals.push_back(0.0);

			const size_t connection_type = matrix[idx][6];

			// save the indices of the nodes that are connected by the considered Link
			vector<size_t> indices = {j, k, connection_type};
			Connections.emplace_back(indices);

			// define equilibrium distances
			// polymer-polymer links
			if (connection_type==0){
				vector<double> eq = {0.0, 0.0, static_cast<double>(1)};
				eq_dist.emplace_back(eq);
			}	
			// filler-filler links
			else if (connection_type==1) {
				double djk_sq = 0.0;
				for (size_t i = 0; i < dim; i++)
				{
					const double lbox = box_size[i];
					double xjk = Nodes[j][i] - Nodes[k][i];
					xjk -= round(xjk / lbox) * lbox;
					djk_sq += xjk * xjk;
				}
				const double djk = sqrt(djk_sq);

				vector<double> eq = {djk, R_factor_ff_sq * djk_sq, static_cast<double>(1)};
				eq_dist.emplace_back(eq);
			}
			// filler-polymer links
			else if (connection_type==2) {
				double djk_sq = 0.0;
				for (size_t i = 0; i < dim; i++)
				{
					const double lbox = box_size[i];
					double xjk = Nodes[j][i] - Nodes[k][i];
					xjk -= round(xjk / lbox) * lbox;
					djk_sq += xjk * xjk;
				}
				const double djk = sqrt(djk_sq);

				vector<double> eq = {djk, R_factor_pf_sq * djk_sq, static_cast<double>(1)};
				eq_dist.emplace_back(eq);
			}

			// save information whether a link crosses the boundary
			vector<double> boundary;
			for (size_t i = 0; i < dim; i++)
			{
				const double lbox = box_size[i];
				double dist = Nodes[j][i] - Nodes[k][i];
				const double boundary_info = round(dist / lbox);
				boundary.push_back(boundary_info);
				dist -= boundary_info * lbox;
			}

			Periodic_BC[j].emplace_back(boundary);

			boundary.clear();
		}

		matrix.clear();
		std::cout << "# Reading " + _initial_linkage_file + " successful." << std::endl;
	}

	// this function reads the type of the nodes from a file, i.e. whether they are filler or not
	void read_initial_filler(const std::string &_initial_filler_file)
	{
		std::cout << "# Reading " + _initial_filler_file + " ..." << std::endl;

		std::ifstream file(_initial_filler_file);
		std::string line;
		std::vector<std::vector<double>> matrix;

		if (!file)
		{
			std::cout << "# Cannot find " + _initial_filler_file << "." << std::endl;
			std::cin.ignore();
			exit(0);
		}

		while (std::getline(file, line))
		{
			std::stringstream ss(line);
			std::vector<double> line_data;

			double value;
			while (ss >> value)
			{
				line_data.push_back(value);
			}
			matrix.push_back(line_data);
		}
		file.close();

		// assign the type ''filler'' = 1 or ''no filler'' = 0 to the nodes
		if (matrix.size() != NodesN)
		{
			std::cout << "# Number of nodes given in the file does not fit the given parameters." << std::endl;
			exit(0);
		}

		Filler.resize(NodesN);

		for (size_t j = 0; j < NodesN; j++)
		{
			if (matrix[j][0] == 0)
			{
				Filler[j] = false;
			}
			else
			{
				Filler[j] = true;
			}
		}
		matrix.clear();

		// add the type of the links to the vector ''Connections''
		// obtain the equilibrium distances of filler-filler and polymer-filler connections

		eq_dist.clear();

		for (size_t l = 0; l < Connections.size(); l++)
		{
			const size_t node_j = Connections[l][0];
			const size_t node_k = Connections[l][1];

			// type of the nodes
			const bool filler_j = Filler[node_j];
			const bool filler_k = Filler[node_k];

			if (filler_j == filler_k)
			{
				// polymer-polymer links are type 0
				if (filler_j == false)
				{
					Connections[l].push_back(0);

					vector<double> eq = {0.0, 0.0, static_cast<double>(1)};
					eq_dist.emplace_back(eq);
				}
				// filler-filler links are type 1
				else
				{
					Connections[l].push_back(1);

					double djk_sq = 0.0;
					for (size_t i = 0; i < dim; i++)
					{
						const double lbox = box_size[i];
						double xjk = Nodes[node_j][i] - Nodes[node_k][i];
						xjk -= round(xjk / lbox) * lbox;
						djk_sq += xjk * xjk;
					}
					const double djk = sqrt(djk_sq);

					vector<double> eq = {djk, R_factor_ff_sq * djk_sq, static_cast<double>(1)};
					eq_dist.emplace_back(eq);
				}
			}
			// polymer-filler links are type 2
			else
			{
				Connections[l].push_back(2);

				double djk_sq = 0.0;
				for (size_t i = 0; i < dim; i++)
				{
					const double lbox = box_size[i];
					double xjk = Nodes[node_j][i] - Nodes[node_k][i];
					xjk -= round(xjk / lbox) * lbox;
					djk_sq += xjk * xjk;
				}
				const double djk = sqrt(djk_sq);

				vector<double> eq = {djk, R_factor_pf_sq * djk_sq, static_cast<double>(1)};
				eq_dist.emplace_back(eq);
			}
		}
		std::cout << "# Reading " + _initial_filler_file + " successful." << std::endl;
	}

	// this function establishes the links between the nodes by using a neighbor list
	// if the network should have a crack, it is set here
	void nodes_neighbor_list()
	{
		LinksN = 0;

		Neighbors.clear();
		Periodic_BC.clear();

		// check if fillers should be embedded in the network
		if (filled == true)
		{
			// set up the list of nearest neighbors for each node
			set_nodes_neighbor_list_filler();
		}
		// set neighbors of nodes if the network should not contain fillers
		else
		{
			vector<vector<size_t>> list;
			list.clear();
			list.resize(NodesN);

			vector<vector<size_t>> list_backup;
			list_backup.clear();
			list_backup.resize(NodesN);

			// vector<size_t> neighbor_no;

			vector<vector<vector<double>>> periodic_list;
			periodic_list.clear();
			periodic_list.resize(NodesN);

			for (size_t j = 0; j < NodesN; j++)
			{
				const bool j_no_filler = (Filler[j] == 0);

				for (size_t k = j + 1; k < NodesN; k++)
				{
					const bool k_no_filler = (Filler[k] == 0);

					vector<double> boundary;
					// compute distance^2 between node j and node k
					double distsq_jk = 0.0;
					for (size_t i = 0; i < dim; i++)
					{
						const double lbox = box_size[i];
						double dist = Nodes[j][i] - Nodes[k][i];
						const double boundary_info = round(dist / lbox);
						boundary.push_back(boundary_info);
						dist -= boundary_info * lbox;
						distsq_jk += dist * dist;
					}

					// check if node j is located within the neighbor radius of node k
					// if both of the nodes are cross links
					if ((j_no_filler && k_no_filler) && distsq_jk < nodes_rcutsq && j != k)
					{
						list[j].push_back(k);
						periodic_list[j].emplace_back(boundary);

						list_backup[j].push_back(k);
						list_backup[k].push_back(j);
					}
					// if at least one of the nodes is filler
					else if ((j_no_filler == false || k_no_filler == false) && distsq_jk < filler_rcutsq && j != k)
					{
						list[j].push_back(k);
						periodic_list[j].emplace_back(boundary);

						list_backup[j].push_back(k);
						list_backup[k].push_back(j);
					}
					boundary.clear();
				}

				Neighbors.push_back(list[j]);
				Periodic_BC.emplace_back(periodic_list[j]);

				LinksN += list[j].size();
			}

			list.clear();
			list_backup.clear();
			periodic_list.clear();
		}

		if (CrackCentered == true)
		{
			insert_crack_centered();
		}
		else
		{
			insert_crack();
		}

		Links_per_Node = 2.0 * static_cast<double>(LinksN) / static_cast<double>(NodesN);
		cout << "# Average Number of Links per Node is " << Links_per_Node << endl;

		cout << "# The total number of generated links is " << LinksN << endl;
	}

	// this function sets up the neighbor lists for the nodes in the case that there are fillers required
	// the neighbors are determined by the minimum distance
	void set_nodes_neighbor_list_filler()
	{
		vector<vector<size_t>> list;
		list.clear();
		list.resize(NodesN);

		vector<vector<size_t>> list_backup;
		list_backup.clear();
		list_backup.resize(NodesN);

		vector<vector<vector<double>>> periodic_list;
		periodic_list.clear();
		periodic_list.resize(NodesN);

		vector<vector<double>> distances;
		distances.clear();
		distances.resize(NodesN);

		// vector that contains indices of the nodes that are consideres as filler
		vector<size_t> filler_list;

		for (size_t j = 0; j < NodesN; j++)
		{

			// number of nearest neighbors for a node that is not a filler
			size_t neigh_no = 2 * dim;

			// check whether node j should be a filler
			if (randomnumber(0.0, 1.0) < filler_frac)
			{

				filler_list.emplace_back(j);

				neigh_no = filler_bonds;
			}

			// vectors to store infomation of bonds
			vector<vector<double>> boundary_bonds = periodic_list[j];
			vector<size_t> nodes_bonds = list[j];
			vector<double> dist_bonds = distances[j];

			double max_dist;
			if (nodes_bonds.size() > 0)
			{
				max_dist = *max_element(dist_bonds.begin(), dist_bonds.end());
			}
			else
			{
				// maximum distance of a neighbor node k to node j which is initalized as distance larger than all possible values
				max_dist = init_length;
			}

			// find nearest neighbors with k < j
			for (size_t k = 0; k < NodesN; k++)
			{
				bool j_neigh_k = false;
				// if (k < j) {
				for (size_t l = 0; l < list[k].size(); l++)
				{
					if (list[k][l] == j)
					{
						// bondsJ += 1;
						j_neigh_k = true;
					}
				}
				//}

				// if j is not a neighbor of node k, check whether it is close enough to be a neighbor
				if (j_neigh_k == false && k != j)
				{
					vector<double> boundary;
					// compute distance^2 between node j and node k
					double distsq_jk = 0.0;
					for (size_t i = 0; i < dim; i++)
					{
						const double lbox = box_size[i];
						double dist = Nodes[j][i] - Nodes[k][i];
						const double boundary_info = round(dist / lbox);
						boundary.push_back(boundary_info);
						dist -= boundary_info * lbox;
						distsq_jk += dist * dist;
					}

					// if there is not already the maximum number of neighbors assigned to node j, add node k to the neighbor list
					if (nodes_bonds.size() < neigh_no)
					{
						dist_bonds.emplace_back(distsq_jk);
						nodes_bonds.emplace_back(k);
						boundary_bonds.emplace_back(boundary);

						max_dist = *max_element(dist_bonds.begin(), dist_bonds.end());
					}
					// if the maximum number of neighbors is already assigned to node j, find the ones with the smallest distance to node j
					else if (distsq_jk < max_dist && nodes_bonds.size() == neigh_no)
					{
						auto replace_index = max_element(dist_bonds.begin(), dist_bonds.end());
						const size_t index = std::distance(dist_bonds.begin(), replace_index);

						dist_bonds[index] = distsq_jk;
						nodes_bonds[index] = k;
						boundary_bonds[index] = boundary;

						max_dist = *max_element(dist_bonds.begin(), dist_bonds.end());
					}
					boundary.clear();
				}
			}

			// now we have the nodes with the smallest distances to node j
			list[j] = nodes_bonds;
			periodic_list[j] = boundary_bonds;
			distances[j] = dist_bonds;
			for (size_t l = 0; l < nodes_bonds.size(); l++)
			{
				const size_t k = nodes_bonds[l];

				list[k].push_back(j);
				periodic_list[k].emplace_back(boundary_bonds[l]);
				distances[k].emplace_back(dist_bonds[l]);
			}
			// cout << list[j].size() << endl;
			nodes_bonds.clear();
			dist_bonds.clear();
			boundary_bonds.clear();
		}

		// the neighbors have to be assigned to node j such that their indices are always larger than j
		// and such that each connection is counted once
		Neighbors.resize(NodesN);
		Periodic_BC.resize(NodesN);
		for (size_t j = 0; j < NodesN; j++)
		{
			for (size_t l = 0; l < list[j].size(); l++)
			{
				const size_t k = list[j][l];
				if (k > j)
				{
					bool neighbor = false;
					for (size_t m = 0; m < Neighbors[j].size(); m++)
					{
						if (k == Neighbors[j][m])
						{
							neighbor = true;
							m = Neighbors[j].size();
						}
					}
					if (neighbor == false)
					{
						Neighbors[j].emplace_back(k);
						Periodic_BC[j].emplace_back(periodic_list[j][l]);
						LinksN += 1;
					}
				}
			}
		}
		// cout << LinksN << endl;
		list.clear();
		list_backup.clear();
		periodic_list.clear();
		distances.clear();
	}

	// morphology generator: assign a fraction filler_frac of nodes to the type ''filler''
	void set_filler_nodes()
	{
		Filler.clear();
		for (size_t j = 0; j < NodesN; j++)
		{
			if (randomnumber(0.0, 1.0) < filler_frac)
			{
				Filler.emplace_back(true);
			}
			else
			{
				Filler.emplace_back(false);
			};
		}
	}

	// compute the total surface energy of the network based on the surface tensions and the network confuguration
	double morphology_energy()
	{
		size_t n_polymer_filler = 0;
		for (size_t j = 0; j < NodesN; j++)
		{
			for (size_t m = 0; m < filler_neighbors[j].size(); m++)
			{									   // Neighbors[j].size(); m++) {
				size_t k = filler_neighbors[j][m]; // Neighbors[j][m];
				if (Filler[j] != Filler[k])
				{
					n_polymer_filler += 1;
				}
			}
		}
		const double W = interface_polymer_filler_factor * static_cast<double>(n_polymer_filler);
		return W;
	}

	// set up a neighbor list for the nodes which is used to generate the filler morphology
	void set_morphology_neighbors()
	{
		filler_neighbors.clear();
		filler_neighbors.resize(NodesN);
		for (size_t j = 0; j < NodesN; j++)
		{
			for (size_t k = j + 1; k < NodesN; k++)
			{
				double rjk_sq = 0.0;
				for (size_t i = 0; i < dim; i++)
				{
					double xjk = Nodes[j][i] - Nodes[k][i];
					const double lbox = box_size[i];
					const double boundary_condition = round(xjk / lbox);
					xjk -= boundary_condition * lbox;
					rjk_sq += xjk * xjk;
				}
				if (rjk_sq < r_morphology_sq)
				{
					filler_neighbors[j].push_back(k);
					filler_neighbors[k].push_back(j);
				}
			}
		}
		std::cout << "# Neighbor list for morhology generator set up." << std::endl;
	}

	void run_morphology_generator()
	{

		set_filler_nodes();

		set_morphology_neighbors();

		/* // save initial filler morphology
		fstream filler_nodes;
		filler_nodes.open("fillernodesInitial.txt", ios::app);
		filler_nodes << "# Number of nodes " << NodesN << endl;
		filler_nodes << "# box length " << box_size[0] << endl;
		filler_nodes << "# surface tension polymer dispersive " << surface_tension_polymer_d << endl;
		filler_nodes << "# surface tension polymer polar " << surface_tension_polymer_p << endl;
		filler_nodes << "# surface tension filler dispersive " << surface_tension_filler_d << endl;
		filler_nodes << "# surface tension filler polar " << surface_tension_filler_p << endl;
		filler_nodes << "# filler?	x	y" << endl;
		for (size_t j = 0; j < NodesN; j++)
		{
			filler_nodes << Filler[j] << "	" << Nodes[j][0] << "	" << Nodes[j][1] << endl;
		}
		filler_nodes.close(); */

		/* size_t step_W_total = 0;
		double W_total = morphology_energy();
		fstream filler_energies;
		filler_energies.open("fillerenergies.txt", ios::app);
		filler_energies << "# MCstep		W_total" << endl;
		filler_energies << step_W_total << "	" << setprecision(10) << W_total << endl;
		step_W_total += 1; */
		mtx.lock();
		// initialize random numbers
		unsigned long long tsc = rdtsc();
		long seed = (long)(tsc & 0xFFFFFFFF);
		idum = -seed;

		// size_t step_config = 5000;

		// MC flocculation algorithm
		for (size_t MCstep = 0; MCstep < MCsteps; MCstep++)
		{

			// pick node at random
			const size_t random_node = round(static_cast<double>(NodesN - 1) * ran2(&idum));
			// select neighbor of this node randomly
			// Note that list_backup contains all pairs of neighbors twice, while Neighbors contains them only once.
			const size_t random_neighbor_idx = round(static_cast<double>(filler_neighbors[random_node].size() - 1) * ran2(&idum));
			const size_t random_neighbor = filler_neighbors[random_node][random_neighbor_idx];

			// check whether these nodes are of different type, otherwise interchanging them would not change a thing
			if (Filler[random_node] != Filler[random_neighbor])
			{
				// identify polymer-filler interfaces within the neighbor lists of the picked nodes
				// Note that it does not matter if an interface is counted twice, since the only affected one
				// / one which might be counted twice is that between the interchanged nodes
				size_t nr_polymer_filler = 0;
				for (size_t m = 0; m < filler_neighbors[random_node].size(); m++)
				{
					const size_t index = filler_neighbors[random_node][m];
					if (Filler[index] != Filler[random_node])
					{
						nr_polymer_filler += 1;
					}
				}
				for (size_t m = 0; m < filler_neighbors[random_neighbor].size(); m++)
				{
					const size_t index = filler_neighbors[random_neighbor][m];
					if (Filler[index] != Filler[random_neighbor])
					{
						nr_polymer_filler += 1;
					}
				}
				double number_polymer_filler = static_cast<double>(nr_polymer_filler);

				// compute original surface energy
				const double W_old = interface_polymer_filler_factor * number_polymer_filler;

				// back up the neighbor lists of the nodes for the interchange
				const vector<size_t> neighbors_node_old = filler_neighbors[random_node];
				const vector<size_t> neighbors_neighbor_old = filler_neighbors[random_neighbor];

				// interchange the neighbor lists of the nodes
				filler_neighbors[random_node] = neighbors_neighbor_old;
				filler_neighbors[random_neighbor] = neighbors_node_old;

				// identify new polymer-filler interfaces and count them
				nr_polymer_filler = 0;
				for (size_t m = 0; m < filler_neighbors[random_node].size(); m++)
				{
					size_t index = filler_neighbors[random_node][m];
					if (index == random_node)
					{
						filler_neighbors[random_node][m] = random_neighbor;
						index = random_neighbor;
					}
					if (Filler[index] != Filler[random_node])
					{
						nr_polymer_filler += 1;
					}
				}
				for (size_t m = 0; m < filler_neighbors[random_neighbor].size(); m++)
				{
					size_t index = filler_neighbors[random_neighbor][m];
					if (index == random_neighbor)
					{
						filler_neighbors[random_neighbor][m] = random_node;
						index = random_node;
					}
					if (Filler[index] != Filler[random_neighbor])
					{
						nr_polymer_filler += 1;
					}
				}
				number_polymer_filler = static_cast<double>(nr_polymer_filler);

				// compute new surface energy
				const double W_new = interface_polymer_filler_factor * number_polymer_filler;

				// evaluate the Metropolis criterion
				if (exp(W_old - W_new) < ran2(&idum))
				{
					// exponential(10, W_old - W_new) < ran2(&idum)){
					//  configuration is rejected

					// reset the interchange of the neighbor lists
					filler_neighbors[random_node] = neighbors_node_old;
					filler_neighbors[random_neighbor] = neighbors_neighbor_old;
				}
				else
				{
					// configuration is accepted

					// interchange the positions
					const vector<double> position_node_old = Nodes[random_node];
					const vector<double> position_neighbor_old = Nodes[random_neighbor];
					Nodes[random_node] = position_neighbor_old;
					Nodes[random_neighbor] = position_node_old;

					// interchange indices of the selected nodes on the neighbor lists of their neighbors
					for (size_t m = 0; m < filler_neighbors[random_node].size(); m++)
					{
						size_t neighbor_m = filler_neighbors[random_node][m];
						if (neighbor_m != random_neighbor)
						{
							for (size_t l = 0; l < filler_neighbors[neighbor_m].size(); l++)
							{
								if (filler_neighbors[neighbor_m][l] == random_neighbor)
								{
									filler_neighbors[neighbor_m][l] = random_node;
									l = filler_neighbors[neighbor_m].size();
								}
							}
						}
					}
					for (size_t m = 0; m < filler_neighbors[random_neighbor].size(); m++)
					{
						size_t neighbor_m = filler_neighbors[random_neighbor][m];
						if (neighbor_m != random_node)
						{
							for (size_t l = 0; l < filler_neighbors[neighbor_m].size(); l++)
							{
								if (filler_neighbors[neighbor_m][l] == random_node)
								{
									filler_neighbors[neighbor_m][l] = random_neighbor;
									l = filler_neighbors[neighbor_m].size();
								}
							}
						}
					}
				}
			}

			/* if (MCstep-1 == step_W_total) {
			   W_total = morphology_energy();
			   filler_energies << MCstep+1 << "	" << setprecision(10) << W_total << endl;
			   if (MCstep < 50) {
				   step_W_total += 1;
			   }
			   else {
				   step_W_total += 50;
			   }
			}  */

			/* if(MCstep == step_config){
				filler_nodes.open("fillernodes"+std::to_string(step_config)+".txt", ios::app);
				filler_nodes << "# Number of nodes " << NodesN << endl;
				filler_nodes << "# box length " << box_size[0] << endl;
				filler_nodes << "# filler fraction " << filler_frac << endl;
				filler_nodes << "# MC steps " << MCsteps << endl;
				filler_nodes << "# surface tension polymer dispersive " << surface_tension_polymer_d << endl;
				filler_nodes << "# surface tension polymer polar " << surface_tension_polymer_p << endl;
				filler_nodes << "# surface tension filler dispersive " << surface_tension_filler_d << endl;
				filler_nodes << "# surface tension filler polar " << surface_tension_filler_p << endl;
				filler_nodes << "# filler?	x	y" << endl;
				for (size_t j = 0; j < NodesN; j++)
				{
					filler_nodes << Filler[j] << "	" << Nodes[j][0] << "	" << Nodes[j][1] << endl;
				}
				filler_nodes.close();

				step_config += 5000;	//10000;	//
			} */
		}
		/* W_total = morphology_energy();
		filler_energies << MCsteps << "	" << setprecision(10) << W_total << endl;

		filler_energies.close(); */

		std::cout << "# Filler morphology generated." << std::endl;
		mtx.unlock();
		// save final filler morphology
		fstream filler_nodes;
		filler_nodes.open("fillernodes.txt", ios::app);
		filler_nodes << "# Number of nodes " << NodesN << endl;
		filler_nodes << "# box length " << box_size[0] << endl;
		filler_nodes << "# filler fraction " << filler_frac << endl;
		filler_nodes << "# MC steps " << MCsteps << endl;
		filler_nodes << "# surface tension polymer dispersive " << surface_tension_polymer_d << endl;
		filler_nodes << "# surface tension polymer polar " << surface_tension_polymer_p << endl;
		filler_nodes << "# surface tension filler dispersive " << surface_tension_filler_d << endl;
		filler_nodes << "# surface tension filler polar " << surface_tension_filler_p << endl;
		filler_nodes << "# filler?	x	y" << endl;
		for (size_t j = 0; j < NodesN; j++)
		{
			filler_nodes << Filler[j] << "	" << Nodes[j][0] << "	" << Nodes[j][1] << endl;
		}
		filler_nodes.close();

		if (dim == 3)
		{
			// save final filler morphology
			filler_nodes.open("fillermorphology.txt", ios::app);
			filler_nodes << "# Number of nodes " << NodesN << endl;
			filler_nodes << "# box length " << box_size[0] << endl;
			filler_nodes << "# filler fraction " << filler_frac << endl;
			filler_nodes << "# MC steps " << MCsteps << endl;
			filler_nodes << "# surface tension polymer dispersive " << surface_tension_polymer_d << endl;
			filler_nodes << "# surface tension polymer polar " << surface_tension_polymer_p << endl;
			filler_nodes << "# surface tension filler dispersive " << surface_tension_filler_d << endl;
			filler_nodes << "# surface tension filler polar " << surface_tension_filler_p << endl;
			filler_nodes << "# filler?	x	y	z" << endl;
			for (size_t j = 0; j < NodesN; j++)
			{
				filler_nodes << Filler[j] << "	" << Nodes[j][0] << "	" << Nodes[j][1] << "	" << Nodes[j][2] << endl;
			}
			filler_nodes.close();
		}

		filler_neighbors.clear();
	}

	// this function inserts a crack into the elastomer network
	void insert_crack()
	{
		CrackArea.clear();
		CrackArea.resize(NodesN);
		std::fill(CrackArea.begin(), CrackArea.end(), false);
		if (CrackCheck == true)
		{

			for (size_t j = 0; j < NodesN; j++)
			{
				for (size_t m = 0; m < Neighbors[j].size(); m++)
				{
					const size_t k = Neighbors[j][m];

					// check if nodes j  and k connected by link l are in the area that is assumed for the crack
					//
					// 	|		|		|	  |		|
					//	o-------o-------o-----o-----o---
					//	|		|		|	  |		|
					//	o-------o-------o-----o-----o---
					//					|	  |		|
					//	o-------o-------o-----o-----o---
					//  |		|		|	  |		|
					//	o-------o-------o-----o-----o---
					// 	|		|		|	  |		|
					//
					// first in x-direction: the crack will be placed around zero along this axis

					const double crack_upper_bound = 0.5 * d_init + 1.0;
					const double crack_lower_bound = -crack_upper_bound;
					const double crack_side_bound = CrackSize * d_init - 0.5 * init_length;

					if (Nodes[j][0] < crack_upper_bound && Nodes[j][0] > crack_lower_bound && Nodes[k][0] > crack_lower_bound && Nodes[k][0] < 0.5 * crack_upper_bound && (Nodes[j][0] < Nodes[k][0] - 0.5 * d_init || Nodes[k][0] < Nodes[j][0] - 0.5 * d_init))
					{
						// second in y- and z-direction: the crack will be placed at the boundary
						if (Nodes[j][1] < crack_side_bound || Nodes[k][1] < crack_side_bound)
						{
							Neighbors[j].erase((Neighbors[j].begin() + m));
							Periodic_BC[j].erase(Periodic_BC[j].begin() + m);
							LinksN -= 1;
							m -= 1;
						}
					}
				}
			}
			std::cout << "# Network got a crack..." << endl;
		}
	}

	// this function inserts a crack in the center of the elastomer network
	void insert_crack_centered()
	{
		const double crack_upper_bound = 0.5 * d_init + 1.0;
		const double crack_lower_bound = -crack_upper_bound;
		const double crack_side_bound_pos = 0.5 * (CrackSize * d_init);
		const double crack_side_bound_neg = -crack_side_bound_pos;

		// the crack will be placed around the origin of the coordinate system
		//
		// 	|		|		|	  |		|		|
		//	o-------o-------o-----o-----o-------o---
		//	|		|		|	  |		|		|
		//	o-------o-------o-----o-----o-------o---
		//	|		|			  		|		|
		//	o-------o-------o-----o-----o-------o---
		//  |		|		|	  |		|		|
		//	o-------o-------o-----o-----o-------o---
		// 	|		|		|	  |		|		|
		//

		CrackArea.clear();
		CrackArea.resize(NodesN);
		std::fill(CrackArea.begin(), CrackArea.end(), false);

		const double xArea = 5.0 * d_init;
		const double yArea = crack_side_bound_pos + 7.0 * d_init;

		for (size_t j = 0; j < NodesN; j++)
		{
			const double xj = Nodes[j][0];
			const double yj = Nodes[j][1];

			for (size_t m = 0; m < Neighbors[j].size(); m++)
			{
				const size_t k = Neighbors[j][m];

				// check if nodes j  and k connected by link l are in the area that is assumed for the crack
				const double xk = Nodes[k][0];
				if (xj < crack_upper_bound && xj > crack_lower_bound && xk > crack_lower_bound && xk < 0.5 * crack_upper_bound && (xj < xk - 0.5 * d_init || xk < xj - 0.5 * d_init))
				{
					const double yk = Nodes[k][1];
					if ((yj < crack_side_bound_pos && yj > crack_side_bound_neg) || (yk < crack_side_bound_pos && yk > crack_side_bound_neg))
					{
						Neighbors[j].erase((Neighbors[j].begin() + m));
						Periodic_BC[j].erase(Periodic_BC[j].begin() + m);
						LinksN -= 1;
						m -= 1;
					}
				}
			}

			// set neighborhood of the crack
			if (xj < xArea && xj > -xArea && yj < yArea && yj > -yArea)
			{
				CrackArea[j] = true;
			}
		}

		/*vector<size_t> CrackNeighbors;
		for (size_t j = 0; j < NodesN; j++) {
			if (CrackArea[j] == true) {
				for (size_t m = 0; m < Neighbors[j].size(); m++) {
					const size_t k = Neighbors[j][m];
					CrackNeighbors.push_back(k);
					for (size_t n = 0; n < Neighbors[k].size(); n++) {
						const size_t l = Neighbors[k][n];
						CrackNeighbors.push_back(l);
						for (size_t o = 0; o < Neighbors[l].size(); o++) {
							const size_t p=Neighbors[l][o];
							CrackNeighbors.push_back(p);
						}
					}
				}
			}
		}

		for (size_t neigh = 0; neigh < CrackNeighbors.size(); neigh++) {
			CrackArea[neigh] = true;
		}*/

		cout << "# Network got a crack in the center ..." << endl;
	}

	// this function sets the number of segments in a link
	// and the defects, i.e. determines which links are  non-crystallizable,
	// and computes the positions of the link centers
	void set_links()
	{
		Segments.clear();
		Segments.resize(NodesN);

		Defects.clear();
		Defects.resize(NodesN);

		Connections.clear();

		Crystals.clear();

		Link_Centers.clear();

		const double mu = d_init * d_init;
		const double sigma = mu / 5.0;

		fstream segments_file;
		segments_file.open("initiallinkageSegments.txt", ios::app);
		// "# link j		neighbor link k		segments n		defect?"

		for (size_t j = 0; j < NodesN; j++)
		{
			const size_t linksJ = Neighbors[j].size();
			for (size_t l = 0; l < linksJ; l++)
			{
				// draw number of segments in a links from Gaussian distribution with mu=<n>=l^2, sigma=mu/5
				double n = round(rnd_gaussian(mu, sigma));
				while (n <= 0.0)
				{
					n = round(rnd_gaussian(mu, sigma));
				}
				Segments[j].push_back(n);

				size_t defect_link;
				// set defects
				if (crystallizing == true && randomnumber(0.0, 1.0) < defect_frac)
				{
					Defects[j].push_back(true);
					defect_link = 1;
				}
				else
				{
					Defects[j].push_back(false);
					defect_link = 0;
				}

				const size_t k = Neighbors[j][l];
				if (Save_Segments == true)
				{
					segments_file << j << "		" << k << "		" << n << "		" << defect_link << endl;
				}

				// set the number of crystals in a link to zero for initial configuration
				Crystals.push_back(0.0);

				// save the indices of the nodes that are connected by the considered Link
				vector<size_t> indices = {j, k};
				Connections.emplace_back(indices);
			}
			// cout << endl;
		}
		segments_file.close();
	}

	// this function sets the number of segments in a link
	// and the defects, i.e. determines which links are  non-crystallizable,
	// and computes the positions of the link centers
	// !!! filler-filler and polymer-filler links are not set up here!!!
	// all links are defined as polmer initially
	void set_initial_links_morphology()
	{
		Segments.clear();
		Segments.resize(NodesN);

		Defects.clear();
		Defects.resize(NodesN);

		Connections.clear();

		Crystals.clear();

		eq_dist.clear();

		const double mu = d_init * d_init;
		const double sigma = mu / 5.0;

		fstream segments_file;
		segments_file.open("initiallinkageSegments.txt", ios::app);
		// "# link j		neighbor link k		segments n		defect?"

		no_ff = 0;
		no_pf = 0;
		for (size_t j = 0; j < NodesN; j++)
		{
			const size_t linksJ = Neighbors[j].size();

			for (size_t l = 0; l < linksJ; l++)
			{
				const size_t k = Neighbors[j][l];

				// draw number of segments in a links from Gaussian distribution with mu=<n>=l^2, sigma=mu/5
				double n = round(rnd_gaussian(mu, sigma));
				while (n <= 0.0)
				{
					n = round(rnd_gaussian(mu, sigma));
				}
				Segments[j].push_back(n);

				if (Filler[j] == Filler[k])
				{
					// both nodes are cross links -> polymer-polymer link
					if (Filler[j] == false)
					{
						// type of link 0=polymer-polymer
						const size_t type = 0;

						size_t defect_link;
						// set defects
						if (crystallizing == true && randomnumber(0.0, 1.0) < defect_frac)
						{
							Defects[j].push_back(true);
							defect_link = 1;
						}
						else
						{
							Defects[j].push_back(false);
							defect_link = 0;
						}

						if (Save_Segments == true)
						{
							segments_file << j << "		" << k << "		" << n << "		" << defect_link << "	" << Filler[j] << "	" << Filler[k] << "	" << type << endl;
						}

						// set the number of crystals in a link to zero for initial configuration
						Crystals.push_back(0.0);

						// save the indices of the nodes that are connected by the considered Link
						// and the type of the link
						vector<size_t> indices = {j, k, type};
						Connections.emplace_back(indices);

						vector<double> eq = {0.0, 0.0, static_cast<double>(1)};
						eq_dist.emplace_back(eq);
					}
					// both of the nodes are filler -> filler-filler link
					else
					{
						// type of the link 1=filler-filler
						const size_t type = 1;

						vector<double> eq = {0.0, 0.0, static_cast<double>(1)};
						eq_dist.emplace_back(eq);

						Defects[j].push_back(true);

						// the filler-filler links are non-crystallizable
						// this is done for computational reasons
						Crystals.push_back(0.0);

						// save the indices of the nodes that are connected by the considered Link
						vector<size_t> indices = {j, k, type};
						Connections.emplace_back(indices);

						if (Save_Segments == true)
						{
							segments_file << j << "		" << k << "		" << 0 << "		" << true << "	" << Filler[j] << "	" << Filler[k] << "	" << type << endl;
						}

						no_ff += 1;
					}
				}
				// the nodes are of different type -> polymer-filler link
				else
				{
					// type of the link 2=polymer-filler
					const size_t type = 2;

					vector<double> eq = {0.0, 0.0, static_cast<double>(1)};
					eq_dist.emplace_back(eq);

					Defects[j].push_back(true);
					// the filler-filler links are non-crystallizable
					// this is done for computational reasons
					Crystals.push_back(0.0);

					// save the indices of the nodes that are connected by the considered Link
					vector<size_t> indices = {j, k, type};
					Connections.emplace_back(indices);

					if (Save_Segments == true)
					{
						segments_file << j << "		" << k << "		" << 0 << "		" << true << "	" << Filler[j] << "	" << Filler[k] << "	" << type << endl;
					}

					no_pf += 1;
				}
			}
		}
		segments_file.close();
	}

	// polymer-polymer links are preserved, but filler-filler and polymer-filler links are modified
	// number of segments is set to 0 and equilibrium distances are defined
	void set_filler_links()
	{
		size_t link_idx = 0;
		for (size_t j = 0; j < NodesN; j++)
		{
			const size_t linksJ = Neighbors[j].size();

			for (size_t l = 0; l < linksJ; l++)
			{
				const size_t k = Neighbors[j][l];

				// if link is of type filler-filler
				if (Filler[j] == Filler[k] && Filler[j] == true)
				{
					// compute current end-to-end distance
					double djk_sq = 0.0;
					for (size_t i = 0; i < dim; i++)
					{
						const double lbox = box_size[i];
						double xjk = Nodes[j][i] - Nodes[k][i];
						xjk -= round(xjk / lbox) * lbox;
						djk_sq += xjk * xjk;
					}

					// compute equilibrium distance based on current end-to-end distance
					const double n_inv = 1.0 / Segments[j][l];
					const double eq_factor_sq = (1.0 - spring_ff_inv * n_inv) * (1.0 - spring_ff_inv * n_inv);
					djk_sq *= eq_factor_sq;

					const double djk = sqrt(djk_sq);

					vector<double> eq = {djk, R_factor_ff_sq * djk_sq, static_cast<double>(1)};
					eq_dist[link_idx] = eq;

					Segments[j][l] = 0.0;
				}

				// if the nodes are of different type -> polymer-filler link
				else if (Filler[j] != Filler[k])
				{
					// compute current end-to-end distance
					double djk_sq = 0.0;
					for (size_t i = 0; i < dim; i++)
					{
						const double lbox = box_size[i];
						double xjk = Nodes[j][i] - Nodes[k][i];
						xjk -= round(xjk / lbox) * lbox;
						djk_sq += xjk * xjk;
					}

					// compute equilibrium distance based on current end-to-end distance
					const double n_inv = 1.0 / Segments[j][l];
					const double eq_factor_sq = (1.0 - spring_pf_inv * n_inv) * (1.0 - spring_pf_inv * n_inv);
					djk_sq *= eq_factor_sq;

					const double djk = sqrt(djk_sq);

					vector<double> eq = {djk, R_factor_pf_sq * djk_sq, static_cast<double>(1)};
					eq_dist[link_idx] = eq;

					Segments[j][l] = 0.0;
				}
				link_idx += 1;
			}
		}
	}

	// this function computes the positions of the centers of the links
	void link_center_positions_morphology()
	{
		//vector<vector<double>> Link_Centers_OLD = Link_Centers;
		size_t m = 0;

		Link_Centers.clear();
		for (size_t j = 0; j < NodesN; j++)
		{
			const size_t linksJ = Neighbors[j].size();
			for (size_t l = 0; l < linksJ; l++)
			{
				const size_t k = Neighbors[j][l];
				vector<double> rjk;

				for (size_t i = 0; i < dim; i++)
				{
					double dr = 0.5 * (Nodes[j][i] + Nodes[k][i]);
					dr -= Periodic_BC[j][l][i] * box_size[i];
					rjk.push_back(dr);
				}
				Link_Centers.emplace_back(rjk);
				m += 1;
			}
		}

		//Link_Centers_OLD.clear();
	}

	// this functions sets up a fixed neighbor list for the link centers
	// it is based on the configuration at the maximum strain
	// in this function different link types are considered
	// only polymer-polymer links are considered
	void links_fixed_neighbor_list_morphology()
	{
		Link_Neighbors.resize(LinksN);

		// compute positions of the nodes at maximum strain, crystals are neglected here
		vector<vector<double>> Nodes_MaxStrain(NodesN, vector<double>(dim));

		vector<double> box_size_Max(dim);
		box_size_Max[0] = init_length * max_strain;
		const double lambda_f_inv = 1.0 / (pow(max_strain, 1.0 / static_cast<double>(dim - 1)));
		for (size_t i = 1; i < dim; i++)
		{
			box_size_Max[i] = init_length * lambda_f_inv;
		}

		for (size_t j = 0; j < NodesN; j++)
		{
			Nodes_MaxStrain[j][0] = Nodes[j][0] * max_strain;
			for (size_t i = 1; i < dim; i++)
			{
				Nodes_MaxStrain[j][i] = Nodes[j][i] * lambda_f_inv;
			}
		}

		// compute positions of link centers at maximum strain based on node positions at maximum strain
		vector<vector<double>> Link_Centers_Max;

		size_t m = 0;
		for (size_t j = 0; j < NodesN; j++)
		{
			const size_t linksJ = Neighbors[j].size();
			for (size_t l = 0; l < linksJ; l++)
			{
				const size_t k = Neighbors[j][l];
				vector<double> rjk;

				for (size_t i = 0; i < dim; i++)
				{
					double dr = 0.5 * (Nodes_MaxStrain[j][i] + Nodes_MaxStrain[k][i]);
					dr -= Periodic_BC[j][l][i] * box_size_Max[i];

					rjk.push_back(dr);
				}

				Link_Centers_Max.emplace_back(rjk);

				m += 1;
			}
		}

		for (size_t l = 0; l < LinksN; l++)
		{
			if (static_cast<size_t>(Connections[l][2]) == 0)
			{
				for (size_t m = l + 1; m < LinksN; m++)
				{
					if (Connections[m][2] == Connections[l][2])
					{
						double rlm_sq = 0.0;
						for (size_t i = 0; i < dim; i++)
						{
							const double lbox = box_size_Max[i];
							double rlm = Link_Centers_Max[l][i] - Link_Centers_Max[m][i];
							rlm -= round(rlm / lbox) * lbox;
							rlm_sq += rlm * rlm;
						}
						if (rlm_sq < rmsq)
						{
							Link_Neighbors[l].push_back(m);
						}
					}
				}
			}
		}
		Link_Centers_Max.clear();
		cout << "# Neighbor list of links is set up." << endl;
	}

	// this function computes force and energy based on the initial configuration of the network,
	// when all links are polymer-polymer
	// !crystal interaction and rupture are ommitted here since the initial configuration is considered!
	void initial_force_and_energy()
	{
		G = 0.0;

		Force.clear();
		Force.resize(NodesN, vector<double>(dim));
		std::for_each(Force.begin(), Force.end(), [](auto &sub)
					  { std::fill(sub.begin(), sub.end(), 0.0); });

		Force_abs.clear();
		Force_abs.resize(NodesN, vector<double>(dim));
		std::for_each(Force_abs.begin(), Force_abs.end(), [](auto &sub)
					  { std::fill(sub.begin(), sub.end(), 0.0); });

		Free_Energy.clear();
		Free_Energy.resize(LinksN);
		std::fill(Free_Energy.begin(), Free_Energy.end(), 0.0);

		// compute free energy contributions of the links themself and the force on the nodes caused by the links
		size_t m = 0;
		for (size_t j = 0; j < NodesN; j++)
		{
			for (size_t l = 0; l < Neighbors[j].size(); l++)
			{
				const size_t k = Neighbors[j][l];
				const double n = Segments[j][l];
				// const double nc = Crystals[m] * lc;

				// vector connecting node j and node k
				vector<double> vec_jk;
				// distance between node j and node k squared
				double distsq_jk = 0.0;

				for (size_t i = 0; i < dim; i++)
				{
					const double coord_j = Nodes[j][i];
					const double coord_k = Nodes[k][i];
					const double boundary_condition = box_size[i] * Periodic_BC[j][l][i];

					double coord_jk_diff = coord_j - coord_k;
					coord_jk_diff -= boundary_condition;

					vec_jk.push_back(coord_jk_diff);
					distsq_jk += coord_jk_diff * coord_jk_diff;
				}

				const double dist_jk = sqrt(distsq_jk);

				const double dist_jk_inv = 1.0 / dist_jk;
				const double dist_jk_diff = dist_jk;			   // - nc;
				const double pre_factor = -3.0 * dist_jk_diff / n; //(n - nc);

				// compute force contribution to force on node j by link m
				//double fsq = 0.0;
				for (size_t i = 0; i < dim; i++)
				{
					const double f = pre_factor * vec_jk[i] * dist_jk_inv;
					Force[j][i] += f;
					Force[k][i] -= f;

					const double bc = Periodic_BC[j][l][i];
					if (bc > 0.0)
					{
						Force_abs[j][i] += f * bc;
					}
					//fsq += f * f;
				}
				
				// compute contribution of conformational entropy to free energy
				const double g = -0.5 * pre_factor * dist_jk_diff;
				Free_Energy[m] += g;
				G += g;
			}
			m += 1;
		}
	}

	// this function computes the forces and free energies
	void force_and_energy(const bool allow_rupture, const bool save_forces, const size_t force_hist_index)
	{
		G = 0.0;

		fstream energy_hist;
		if (save_forces == true)
		{
			energy_hist.open("energyhist" + std::to_string(force_hist_index) + ".txt", ios::app);
			energy_hist << "# number of Links = " << LinksN << endl;
			energy_hist << "# labmda = " << box_size[0] / init_length << endl;
		}
		/* fstream force_hist_area;
		if (save_forces == true)
		{
			force_hist_area.open("forcehistarea" + std::to_string(force_hist_index) + ".txt", ios::app);
			force_hist_area << "# number of Links = " << LinksN << endl;
			force_hist_area << "# labmda = " << box_size[0] / init_length << endl;
		} */

		fstream n_hist;
		if (save_forces == true)
		{
			n_hist.open("nhist" + std::to_string(force_hist_index) + ".txt", ios::app);
			n_hist << "# number of Links = " << LinksN << endl;
			n_hist << "# labmda = " << box_size[0] / init_length << endl;
		}

		fstream r_hist;
		if (save_forces == true)
		{
			r_hist.open("rhist" + std::to_string(force_hist_index) + ".txt", ios::app);
			r_hist << "# number of Links = " << LinksN << endl;
			r_hist << "# labmda = " << box_size[0] / init_length << endl;
		}

		Force.clear();
		Force.resize(NodesN, vector<double>(dim));
		std::for_each(Force.begin(), Force.end(), [](auto &sub)
					  { std::fill(sub.begin(), sub.end(), 0.0); });

		Force_abs.clear();
		Force_abs.resize(NodesN, vector<double>(dim));
		std::for_each(Force_abs.begin(), Force_abs.end(), [](auto &sub)
					  { std::fill(sub.begin(), sub.end(), 0.0); });

		Free_Energy.clear();
		Free_Energy.resize(LinksN);
		std::fill(Free_Energy.begin(), Free_Energy.end(), 0.0);

		// compute contribution of interaction of the semi-crystalline links to the free energy
		for (size_t m = 0; m < LinksN; m++)
		{
			// only polymer-polymer links (type 0) are able to crystallize
			// therefore, only the crystalline segments in them can interact
			if (static_cast<size_t>(Connections[m][2]) == 0)
			{
				const double c_m = Crystals[m];
				if (c_m >= 0.0)
				{
					const double first_term = c_m * theta0_c;
					const double pre_factor = c_m * eta_c;
					Free_Energy[m] += first_term;
					for (size_t k = 0; k < Link_Neighbors[m].size(); k++)
					{
						const size_t l = Link_Neighbors[m][k];
						const double c_l = Crystals[l];

						double distsq_ml = 0.0;
						for (size_t i = 0; i < dim; i++)
						{
							const double lbox = box_size[i];
							double coord_diff = Link_Centers[m][i] - Link_Centers[l][i];
							coord_diff -= round(coord_diff / lbox) * lbox;
							distsq_ml += coord_diff * coord_diff;
						}

						if (distsq_ml < rcutsq)
						{
							const double dist_ml = sqrt(distsq_ml);
							double dg = pre_factor * c_l * exp(-dist_ml * r0_inv);
							Free_Energy[m] -= dg;
							Free_Energy[l] -= dg;
						}
					}
				}
			}
		}

		// compute free energy contributions of the links themself and the force on the nodes caused by the links
		size_t m = 0;
		LinkCheck.clear();
		for (size_t j = 0; j < NodesN; j++)
		{
			for (size_t l = 0; l < Neighbors[j].size(); l++)
			{
				// polymer-polymer link
				if (static_cast<size_t>(Connections[m][2]) == 0)
				{
					const size_t k = Neighbors[j][l];
					const double n = Segments[j][l];

					const double nc = Crystals[m] * lc;

					// vector connecting node j and node k
					vector<double> vec_jk;
					// distance between node j and node k squared
					double distsq_jk = 0.0;

					for (size_t i = 0; i < dim; i++)
					{
						const double coord_j = Nodes[j][i];
						const double coord_k = Nodes[k][i];
						const double boundary_condition = box_size[i] * Periodic_BC[j][l][i];

						double coord_jk_diff = coord_j - coord_k;
						coord_jk_diff -= boundary_condition;

						vec_jk.push_back(coord_jk_diff);
						distsq_jk += coord_jk_diff * coord_jk_diff;
					}

					const double dist_jk = sqrt(distsq_jk);

					if (dist_jk <= n)
					{
						LinkCheck.emplace_back(false);
					}
					else
					{
						LinkCheck.emplace_back(true);
					}

					// std::cout << dist_jk << "		" << n << endl;

					const double dist_jk_inv = 1.0 / dist_jk;
					const double dist_jk_diff = dist_jk - nc;
					const double pre_factor = -3.0 * dist_jk_diff / (n - nc);

					// compute force contribution to force on node j by link m
					// double fsq = 0.0;
					for (size_t i = 0; i < dim; i++)
					{
						const double f = pre_factor * vec_jk[i] * dist_jk_inv;
						Force[j][i] += f;
						Force[k][i] -= f;

						const double bc = Periodic_BC[j][l][i];
						if (bc > 0.0)
						{
							Force_abs[j][i] += f * bc;
						}
						// fsq += f * f;
					}

					if (save_forces == true)
					{
						/* force_hist << fsq << endl;
						if (CrackArea[j] == true)
						{
							force_hist_area << fsq << endl;
						} */
						const double free_chain_length = n - nc;
						n_hist << free_chain_length << endl;
						const double rn_frac = dist_jk / n;
						r_hist << rn_frac << endl;
					}

					// compute contribution of conformational entropy to free energy
					const double g = -0.5 * pre_factor * dist_jk_diff;
					Free_Energy[m] += g;
					G += g;
				}
				// filler-filler link
				else if (static_cast<size_t>(Connections[m][2]) == 1 && eq_dist[m][2] > 0.0)
				{
					const size_t k = Neighbors[j][l];
					// vector connecting node j and node k
					vector<double> vec_jk;
					// distance between node j and node k squared
					double distsq_jk = 0.0;

					for (size_t i = 0; i < dim; i++)
					{
						const double coord_j = Nodes[j][i];
						const double coord_k = Nodes[k][i];
						const double boundary_condition = box_size[i] * Periodic_BC[j][l][i];

						double coord_jk_diff = coord_j - coord_k;
						coord_jk_diff -= boundary_condition;

						vec_jk.push_back(coord_jk_diff);
						distsq_jk += coord_jk_diff * coord_jk_diff;
					}

					const double dist_jk = sqrt(distsq_jk);
					const double elongation = eq_dist[m][0] - dist_jk;
					if (elongation != 0.0)
					{
						const double pre_factor = spring_ff * elongation * 1.0 / dist_jk;

						// compute force contribution to force on node j by link m
						//double fsq = 0.0;
						for (size_t i = 0; i < dim; i++)
						{
							const double f = pre_factor * vec_jk[i];
							Force[j][i] += f;
							Force[k][i] -= f;

							const double bc = Periodic_BC[j][l][i];

							if (bc > 0.0)
							{
								Force_abs[j][i] += f * bc;
							}
							//fsq += f * f;
						}

						// compute contribution to free energy
						const double g = spring_ff_half * elongation * elongation;
						Free_Energy[m] += g;
						G += g;
					}
				}
				// polymer-filler link
				else //if (eq_dist[m][2] > 0.0)
				{
					const size_t k = Neighbors[j][l];
					// vector connecting node j and node k
					vector<double> vec_jk;
					// distance between node j and node k squared
					double distsq_jk = 0.0;

					for (size_t i = 0; i < dim; i++)
					{
						const double coord_j = Nodes[j][i];
						const double coord_k = Nodes[k][i];
						const double boundary_condition = box_size[i] * Periodic_BC[j][l][i];

						double coord_jk_diff = coord_j - coord_k;
						coord_jk_diff -= boundary_condition;

						vec_jk.push_back(coord_jk_diff);
						distsq_jk += coord_jk_diff * coord_jk_diff;
					}

					const double dist_jk = sqrt(distsq_jk);
					const double elongation = eq_dist[m][0] - dist_jk;

					if (elongation != 0.0)
					{
						if (eq_dist[m][2] > 0.0)//(dist_jk < R_factor_pf * eq_dist[m][0])
						{
							const double pre_factor = spring_pf * elongation * 1.0 / dist_jk;

							// compute force contribution to force on node j by link m
							//double fsq = 0.0;
							for (size_t i = 0; i < dim; i++)
							{
								const double f = pre_factor * vec_jk[i];
								Force[j][i] += f;
								Force[k][i] -= f;

								const double bc = Periodic_BC[j][l][i];
								if (bc > 0.0)
								{
									Force_abs[j][i] += f * bc;
								}
								//fsq += f * f;
							}

							// compute contribution to free energy
							const double g = spring_pf_half * elongation * elongation;
							Free_Energy[m] += g;
							G += g;
						}
						else
						{
							const double pre_factor = spring_pf_weak * elongation * 1.0 / dist_jk;

							// compute force contribution to force on node j by link m
							//double fsq = 0.0;
							for (size_t i = 0; i < dim; i++)
							{
								const double f = pre_factor * vec_jk[i];
								Force[j][i] += f;
								Force[k][i] -= f;

								const double bc = Periodic_BC[j][l][i];
								if (bc > 0.0)
								{
									Force_abs[j][i] += f * bc;
								}
								//fsq += f * f;
							}

							// compute contribution to free energy
							const double g = spring_pf_weak_half * elongation * elongation;
							Free_Energy[m] += g;
							G += g;
						}
					}
				}
				m += 1;
			}
		}

		// total free energy of each link is comuted at this point
		// -> rupture criterion based on free energy is checked

		if (allow_rupture == true)
		{
			size_t m = 0;
			for (size_t j = 0; j < NodesN; j++)
			{
				for (size_t l = 0; l < Neighbors[j].size(); l++)
				{
					if (static_cast<size_t>(Connections[m][2]) == 0)
					{
						const double n = Segments[j][l];
						const double c = Crystals[m];
						const double nc = c * lc;

						const double link_crystallinity = nc / n;
						const double g_n = Free_Energy[m] / n;

						// check rupture of the link m
						const bool rupture = link_rupture_morphology(g_n, j, l, m, link_crystallinity);
						if (rupture)
						{
							const size_t k = Neighbors[j][l];
							// distance between node j and node k squared
							double distsq_jk = 0.0;

							for (size_t i = 0; i < dim; i++)
							{
								const double coord_j = Nodes[j][i];
								const double coord_k = Nodes[k][i];
								const double boundary_condition = box_size[i] * Periodic_BC[j][l][i];

								double coord_jk_diff = coord_j - coord_k;
								coord_jk_diff -= boundary_condition;

								distsq_jk += coord_jk_diff * coord_jk_diff;
							}

							const double dist_jk = sqrt(distsq_jk);
							LinksRuptured << n << "			" << dist_jk << "			" << c << "		" << nc << "		" << nc / n << "		" << Defects[j][l] << "		" << dist_jk / eq_dist[m][0] << endl;
							l -= 1;
							m -= 1;
						}
					}
					m += 1;
				}
			}
		}

		if (save_forces == true)
		{
			// free energy histogram is saved

			size_t m = 0;
			for (size_t j = 0; j < NodesN; j++)
			{
				for (size_t l = 0; l < Neighbors[j].size(); l++)
				{
					if (static_cast<size_t>(Connections[m][2]) == 0)
					{
						const double g_n = Free_Energy[m] / Segments[j][l];
						energy_hist << g_n << endl;
					}
					m += 1;
				}
			}

			energy_hist.close();
			// force_hist_area.close();
			n_hist.close();
			r_hist.close();
		}
	}

	// if rupture of links is enabled, this function checks whether a link is broken
	bool link_rupture_morphology(const double _g_n, const size_t node_index, const size_t neighbor_index, const size_t link_index, const double _link_crystallinity)
	{
		bool rupture = false;
		if (RuptureCheck && _link_crystallinity <= RuptureCryst)
		{
			if (_g_n > RuptureEnergy)
			{
				Neighbors[node_index].erase(Neighbors[node_index].begin() + neighbor_index);

				Periodic_BC[node_index].erase(Periodic_BC[node_index].begin() + neighbor_index);
				Segments[node_index].erase(Segments[node_index].begin() + neighbor_index);
				Defects[node_index].erase(Defects[node_index].begin() + neighbor_index);

				Connections.erase(Connections.begin() + link_index);
				Crystals.erase(Crystals.begin() + link_index);
				Free_Energy.erase(Free_Energy.begin() + link_index);

				Link_Centers.erase(Link_Centers.begin() + link_index);

				eq_dist.erase(eq_dist.begin() + link_index);

				for (size_t l = 0; l < LinksN; l++)
				{
					for (size_t m = 0; m < Link_Neighbors[l].size(); m++)
					{
						if (l < link_index && Link_Neighbors[l][m] == link_index)
						{
							Link_Neighbors[l].erase(Link_Neighbors[l].begin() + m);
							m -= 1;
						}
						else if (Link_Neighbors[l][m] > link_index)
						{
							Link_Neighbors[l][m] -= 1;
						}
					}
				}
				Link_Neighbors.erase(Link_Neighbors.begin() + link_index);

				rupture = true;
				LinksN -= 1;
				LinksCracked += 1;
				if (LinkCheck[link_index] == true)
				{
					LinksCrackedRN += 1;
				}

				// cout << "A Link cracked." << endl;
			}
		}
		return rupture;
	}

	// this function minimizes the free energy of the initial system
	// which only contains polymer-polymer links
	// it works analogous to the conventionally applied FIRE algorithm function,
	// but energy and force computation does not consider the type of the chains
	void initial_FIRE_optimizer()
	{
		vector<vector<double>> velocity(NodesN, vector<double>(dim, 0.0));
		vector<vector<double>> acceleration(NodesN, vector<double>(dim, 0.0));

		double alpha = FIRE_alpha_start;
		double dt = FIRE_dt_start;
		double half_dt = 0.5 * dt;
		double dt_inv = 1.0 / dt;
		double alpha_red = 1.0 - alpha;
		size_t it = 0;
		bool convergence = false;

		for (size_t iter = 0; iter < FIRE_N_max; iter++)
		{
			double P = 0.0;
			double G_old = G;

			bool update_links_BC = false;

			// "velocity"-Verlet algorithm to obtain positions and velocities of the nodes
			for (size_t j = 0; j < NodesN; j++)
			{
				for (size_t i = 0; i < dim; i++)
				{
					const double a = acceleration[j][i];

					Nodes[j][i] += dt * velocity[j][i] + half_dt * dt * a;
					const double lbox = box_size[i];

					// check if node crossed boundary of simulation box
					const double boundary_check = round(Nodes[j][i] / lbox);

					// if the node crossed the boundary of the original simulation box, project it back into it
					if ((boundary_check) != 0.0)
					{
						Nodes[j][i] -= boundary_check * lbox;
						update_links_BC = true;
					}
					velocity[j][i] += half_dt * a;
				}
			}

			if (update_links_BC == true)
			{
				for (size_t j = 0; j < NodesN; j++)
				{
					const size_t linksJ = Neighbors[j].size();
					for (size_t l = 0; l < linksJ; l++)
					{
						const size_t k = Neighbors[j][l];
						for (size_t i = 0; i < dim; i++)
						{
							const double dist = Nodes[j][i] - Nodes[k][i];
							const double boundary = round(dist / box_size[i]);
							Periodic_BC[j][l][i] = boundary;
						}
					}
				}
				update_links_BC = false;
			}

			link_center_positions_morphology();

			initial_force_and_energy();

			for (size_t j = 0; j < NodesN; j++)
			{
				double v_abs_sq = 0.0;
				double f_abs_sq = 0.0;
				for (size_t i = 0; i < dim; i++)
				{

					const double f = Force[j][i];
					const double a = f * FIRE_mass_inv;
					acceleration[j][i] = a;
					const double v = velocity[j][i] + half_dt * a;
					velocity[j][i] = v;

					P += f * v;

					v_abs_sq += v * v;
					f_abs_sq += f * f;
				}

				const double v_abs = sqrt(v_abs_sq);
				const double f_abs_inv = 1.0 / sqrt(f_abs_sq);

				for (size_t i = 0; i < dim; i++)
				{
					const double f = Force[j][i];

					double v_FIRE = alpha_red * velocity[j][i];
					if (f != 0.0)
					{
						v_FIRE += alpha * f * f_abs_inv * v_abs;
					}

					velocity[j][i] = v_FIRE;
				}
			}

			if (P < 0.0)
			{
				dt = dt * FIRE_f_dec;
				alpha = FIRE_alpha_start;
				it = iter;

				half_dt = 0.5 * dt;
				dt_inv = 1.0 / dt;
				alpha_red = 1.0 - alpha;

				std::for_each(velocity.begin(), velocity.end(), [](auto &sub)
							  { std::fill(sub.begin(), sub.end(), 0.0); });
			}
			else if (P >= 0.0 && (iter - it) > FIRE_N_min)
			{
				dt = std::min(dt * FIRE_f_inc, FIRE_dt_max);
				alpha *= FIRE_f_alpha;

				half_dt = 0.5 * dt;
				dt_inv = 1.0 / dt;
				alpha_red = 1.0 - alpha;
			}

			const double G_diff = abs(G * 1.0 / G_old - 1.0) * dt_inv;

			if (iter > 1 && G_diff < FIRE_w_rel)
			{
				cout << "# Global energy minimization has finished after " << iter << " iterations..." << endl;
				iter = FIRE_N_max;
				convergence = true;
			}
		}
		if (convergence == false)
		{
			cout << "# FIRE algorithm for global energy minimization does not converge!" << endl;
		}
	}

	// this function minimizes the total free energy of the system
	// by using the "velocity"-Verlet algorithm to obtain the positions and velocities of the nodes
	// and the FIRE algorithm
	void FIRE_optimizer_morphology()
	{
		vector<vector<double>> velocity(NodesN, vector<double>(dim, 0.0));
		vector<vector<double>> acceleration(NodesN, vector<double>(dim, 0.0));

		double alpha = FIRE_alpha_start;
		double dt = FIRE_dt_start;
		double half_dt = 0.5 * dt;
		double dt_inv = 1.0 / dt;
		double alpha_red = 1.0 - alpha;
		size_t it = 0;
		bool convergence = false;

		for (size_t iter = 0; iter < FIRE_N_max; iter++)
		{
			double P = 0.0;
			double G_old = G;

			bool update_links_BC = false;

			// "velocity"-Verlet algorithm to obtain positions and velocities of the nodes
			for (size_t j = 0; j < NodesN; j++)
			{
				for (size_t i = 0; i < dim; i++)
				{
					const double a = acceleration[j][i];

					Nodes[j][i] += dt * velocity[j][i] + half_dt * dt * a;
					const double lbox = box_size[i];

					// check if node crossed boundary of simulation box
					const double boundary_check = round(Nodes[j][i] / lbox);

					// if the node crossed the boundary of the original simulation box, project it back into it
					if ((boundary_check) != 0.0)
					{
						Nodes[j][i] -= boundary_check * lbox;
						update_links_BC = true;
					}
					velocity[j][i] += half_dt * a;
				}
			}

			if (update_links_BC == true)
			{
				for (size_t j = 0; j < NodesN; j++)
				{
					const size_t linksJ = Neighbors[j].size();
					for (size_t l = 0; l < linksJ; l++)
					{
						const size_t k = Neighbors[j][l];
						for (size_t i = 0; i < dim; i++)
						{
							const double dist = Nodes[j][i] - Nodes[k][i];
							const double boundary = round(dist / box_size[i]);
							Periodic_BC[j][l][i] = boundary;
						}
					}
				}
				update_links_BC = false;
			}

			link_center_positions_morphology();

			force_and_energy(false, false, 0);

			for (size_t j = 0; j < NodesN; j++)
			{
				double v_abs_sq = 0.0;
				double f_abs_sq = 0.0;
				for (size_t i = 0; i < dim; i++)
				{

					const double f = Force[j][i];
					const double a = f * FIRE_mass_inv;
					acceleration[j][i] = a;
					const double v = velocity[j][i] + half_dt * a;
					velocity[j][i] = v;

					P += f * v;

					v_abs_sq += v * v;
					f_abs_sq += f * f;
				}

				const double v_abs = sqrt(v_abs_sq);
				const double f_abs_inv = 1.0 / sqrt(f_abs_sq);

				for (size_t i = 0; i < dim; i++)
				{
					const double f = Force[j][i];
					// const double v_FIRE = alpha_red * velocity[j][i] + alpha * f * f_abs_inv * v_abs;
					double v_FIRE = alpha_red * velocity[j][i];
					if (f != 0.0)
					{
						v_FIRE += alpha * f * f_abs_inv * v_abs;
					}

					velocity[j][i] = v_FIRE;
				}
			}

			// cout << P << endl;

			if (P < 0.0)
			{
				dt = dt * FIRE_f_dec;
				alpha = FIRE_alpha_start;
				it = iter;

				half_dt = 0.5 * dt;
				dt_inv = 1.0 / dt;
				alpha_red = 1.0 - alpha;

				std::for_each(velocity.begin(), velocity.end(), [](auto &sub)
							  { std::fill(sub.begin(), sub.end(), 0.0); });
			}
			else if (P >= 0.0 && (iter - it) > FIRE_N_min)
			{
				dt = std::min(dt * FIRE_f_inc, FIRE_dt_max);
				alpha *= FIRE_f_alpha;

				half_dt = 0.5 * dt;
				dt_inv = 1.0 / dt;
				alpha_red = 1.0 - alpha;
			}

			const double G_diff = abs(G * 1.0 / G_old - 1.0) * dt_inv;
			if (iter > 1 && G_diff < FIRE_w_rel)
			{
				cout << "# Global energy minimization has finished after " << iter << " iterations..." << endl;
				iter = FIRE_N_max;
				convergence = true;
			}

		}

		if (convergence == false)
		{
			cout << "# FIRE algorithm for global energy minimization does not converge!" << endl;
		}
	}

	// this function checks whether a filler-filler or a polymer-filler links exceeds the cut-off elongation or not
	// it is called after energy minimization by the FIRE algorithm
	// the number of filler-filler and polymer-filler bonds is counted
	void filler_reversible_bonds()
	{
		no_ff = 0;
		no_pf = 0;
		// check if cut-off of polymer-filler or filler-filler bonds is enabled
		if (cut_ff || cut_pf)
		{
			size_t m = 0;
			for (size_t j = 0; j < Nodes.size(); j++)
			{
				for (size_t l = 0; l < Neighbors[j].size(); l++)
				{
					// if cut-off of filler-filler bonds is enabled and the considered link is a filler-filler bond,
					// check if cut-off distance is exeeded
					if (static_cast<size_t>(Connections[m][2]) == 1)
					{
						if (cut_ff)
						{
							const size_t k = Neighbors[j][l];

							double distjk_sq = 0.0;
							for (size_t i = 0; i < dim; i++)
							{
								double rjk = Nodes[j][i] - Nodes[k][i];
								rjk -= Periodic_BC[j][l][i] * box_size[i];
								distjk_sq += rjk * rjk;
							}

							// if cut-off distance is exceeded, the interaction is cut off
							if (distjk_sq > eq_dist[m][1])
							{
								eq_dist[m][2] = 0.0;
							}
							else
							{
								eq_dist[m][2] = 1.0;
								no_ff += 1;
							}
						}
						else
						{
							no_ff += 1;
						}
					}
					// if cut-off of polymer-filler bonds is enabled and the considered link is a polymer-filler bond,
					// check if cut-off distance is exeeded
					else if (static_cast<size_t>(Connections[m][2]) == 2)
					{
						if (cut_pf)
						{
							const size_t k = Neighbors[j][l];

							double distjk_sq = 0.0;
							for (size_t i = 0; i < dim; i++)
							{
								double rjk = Nodes[j][i] - Nodes[k][i];
								rjk -= Periodic_BC[j][l][i] * box_size[i];
								distjk_sq += rjk * rjk;
							}

							// if cut-off distance is exceeded, the interaction is cut off
							if (distjk_sq > eq_dist[m][1])
							{
								eq_dist[m][2] = 0.0;
							}
							else
							{
								eq_dist[m][2] = 1.0;
								no_pf += 1;
							}
						}
						else
						{
							no_pf += 1;
						}
					}
					m += 1;
				}
			}
		}
		else
		{
			for (size_t m = 0; m < LinksN; m++)
			{
				if (static_cast<size_t>(Connections[m][2]) == 1)
				{
					no_ff += 1;
				}
				else if (static_cast<size_t>(Connections[m][2]) == 2)
				{
					no_pf += 1;
				}
			}
		}
	}

	// this function computes the free energy of a link -> used for crystallization
	// the arguments of this function are the number of crystals in the link,
	// the number of segments in the link, the index of the link
	// and the indices of the nodes that are connected by the link
	double local_free_energy(const double c, const double n, const size_t link_index,
							 const size_t j, const size_t k, const size_t neigh_index)
	{
		const double clc = c * lc;
		double risq = 0.0;
		for (size_t i = 0; i < dim; i++)
		{
			double dr = Nodes[j][i] - Nodes[k][i];
			dr -= box_size[i] * Periodic_BC[j][neigh_index][i];
			risq += dr * dr;
		}
		const double ri_red = sqrt(risq) - clc;
		double g = 1.5 * ri_red * ri_red / (n - clc) + theta0 * clc;

		double interaction = 0.0;
		for (size_t l = 0; l < link_index; l++)
		{
			for (size_t m = 0; m < Link_Neighbors[l].size(); m++)
			{
				const size_t neighbor_link = Link_Neighbors[l][m];
				if (neighbor_link == link_index)
				{
					double distsq = 0.0;
					for (size_t i = 0; i < dim; i++)
					{
						const double lbox = box_size[i];
						double dr = Link_Centers[link_index][i] - Link_Centers[l][i];
						dr -= round(dr / lbox) * lbox;
						distsq += dr * dr;
					}
					if (distsq < rcutsq)
					{
						const double dist = sqrt(distsq);
						interaction += Crystals[l] * exp(-dist * r0_inv);
					}
				}
			}
		}

		for (size_t m = 0; m < Link_Neighbors[link_index].size(); m++)
		{
			const size_t neighbor_index = Link_Neighbors[link_index][m];
			double distsq = 0.0;
			for (size_t i = 0; i < dim; i++)
			{
				const double lbox = box_size[i];
				double dr = Link_Centers[link_index][i] - Link_Centers[neighbor_index][i];
				dr -= round(dr / lbox) * lbox;
				distsq += dr * dr;
			}
			if (distsq < rcutsq)
			{
				const double dist = sqrt(distsq);
				interaction += Crystals[neighbor_index] * exp(-dist * r0_inv);
			}
		}

		interaction *= clc * eta;
		g -= interaction;
		return g;
	}

	// this function checks whether crystallization is preferred by checking local free energy
	void crystallization_check()
	{
		if (crystallizing == true)
		{
			size_t l = 0;
			for (size_t j = 0; j < NodesN; j++)
			{
				for (size_t m = 0; m < Neighbors[j].size(); m++)
				{
					const double g_OLD = Free_Energy[l];
					const double c_OLD = Crystals[l];

					if (c_OLD >= 0 && Defects[j][m] == false)
					{
						const size_t k = Neighbors[j][m];
						const double n = Segments[j][m];
						const double c_p1 = c_OLD + 1.0;
						const double g_p1 = local_free_energy(c_p1, n, l, j, k, m);

						const double dwp = g_p1 - g_OLD;

						if (dwp < 0.0 && (n - c_p1 * lc) >= 1.0)
						{
							Crystals[l] += 1.0;
						}
						else if (c_OLD >= 1.0)
						{
							const double c_m1 = c_OLD - 1.0;
							const double g_m1 = local_free_energy(c_m1, n, l, j, k, m);
							const double dwm = g_m1 - g_OLD;
							if (dwm < 0.0)
							{
								Crystals[l] -= 1.0;
							}
						}
					}
					l += 1;
				}
			}

			cout << "# Creating and melting of crystals finished..." << endl;
		}
	}

	// this function stretches the simulation box
	// and applies a corresponding affine transformation to the positions of the nodes
	void stretching(const double _strain, const double _previous_strain)
	{
		const double lambda_f = _strain;
		// stretch the simulation box along direction 0 by conserving the volume
		box_size[0] = init_length * lambda_f;
		// cout << box_size[0] << endl;
		const double lambda_f_inv = 1.0 / (pow(lambda_f, 1.0 / static_cast<double>(dim - 1)));
		for (size_t i = 1; i < dim; i++)
		{
			box_size[i] = init_length * lambda_f_inv;
			// cout << box_size[i] << endl;
		}

		// apply affine displacements to node positions
		if (lambda_f > 1.0)
		{
			const double scaling_lambda = lambda_f / _previous_strain;
			const double scaling_lambda_inv = 1.0 / (pow(scaling_lambda, 1.0 / static_cast<double>(dim - 1)));
			for (size_t j = 0; j < NodesN; j++)
			{
				Nodes[j][0] *= scaling_lambda;
				// cout << Nodes[j][0] << endl;
				for (size_t i = 1; i < dim; i++)
				{
					Nodes[j][i] *= scaling_lambda_inv;
				}
			}
		}

		// compute A1/A2
		// 2D: A1/A2 = lambda^-2; 3D: A1/A2 = lambda^-3/2
		ratio_A1A2 = 1.0 / lambda_f * lambda_f_inv;
	}

	void stress()
	{
		// double sigma1 = 0.0;
		// double sigma2 = 0.0;
		sigma1 = 0.0;
		sigma2 = 0.0;
		for (size_t j = 0; j < NodesN; j++)
		{
			sigma1 += Force_abs[j][0];
			sigma2 += Force_abs[j][1];
			// cout << Force[j][0] << "		" << Force[j][1] << endl;
		}
		sigma = A0_inv * (sigma1 - ratio_A1A2 * sigma2);
		std::cout << "# sigma = " << sigma << std::endl;
	}

	void crystallinity_index()
	{
		double chi_enum = 0.0;
		double chi_denom = 0.0;

		for (size_t l = 0; l < LinksN; l++)
		{
			chi_enum += Crystals[l];

			const size_t j = Connections[l][0];
			const size_t k = Connections[l][1];

			for (size_t m = 0; m < Neighbors[j].size(); m++)
			{
				if (Neighbors[j][m] == k)
				{
					chi_denom += Segments[j][m];
				}
			}
		}
		chi_enum *= lc;
		chi = chi_enum / chi_denom;
	}

	void save_data()
	{
		if (save == true)
		{
			// save positions of the nodes
			// and size of the simulation box
			savePos.open("pos" + std::to_string(save_counter) + ".txt", ios::app);
			saveBox.open("box" + std::to_string(save_counter) + ".txt", ios::app);
			for (size_t i = 0; i < dim; i++)
			{
				for (size_t j = 0; j < NodesN; j++)
				{
					savePos << Nodes[j][i] << endl;
				}
				saveBox << box_size[i] << endl;
			}
			savePos.close();
			saveBox.close();

			// save information about crystallization of a link
			saveCryst.open("cryst" + std::to_string(save_counter) + ".txt", ios::app);
			saveCryst << "# j	k	c	connected?	cl_c	n" << endl;
			for (size_t l = 0; l < LinksN; l++)
			{
				const size_t j = Connections[l][0];
				const size_t k = Connections[l][1];
				const double cryst = Crystals[l];
				const double link_cryst = cryst * lc;

				double seg;
				for (size_t m = 0; m < Neighbors[j].size(); m++)
				{
					if (Neighbors[j][m] == k)
					{
						seg = Segments[j][m];
						m = Neighbors[j].size();
					}
				}

				saveCryst << j << "		" << k << "		" << cryst << "		" << eq_dist[l][2] << "		" << link_cryst << "		" << seg << endl;
			}
			saveCryst.close();

			if (local_link_data_save == true)
			{
				save_all_local_data(save_counter);
			}

			save_counter += 1;
		}
	}

	void save_initial_r_of_links()
	{
		size_t link_id = 0;
		for (size_t j = 0; j < NodesN; j++)
		{
			const size_t linksJ = Neighbors[j].size();

			for (size_t l = 0; l < linksJ; l++)
			{
				const size_t k = Neighbors[j][l];

				// if link is of type polymer-polymer
				if (Filler[j] == Filler[k] && Filler[j] == false)
				{
					// compute current end-to-end distance
					double djk_sq = 0.0;
					for (size_t i = 0; i < dim; i++)
					{
						const double lbox = box_size[i];
						double xjk = Nodes[j][i] - Nodes[k][i];
						xjk -= round(xjk / lbox) * lbox;
						djk_sq += xjk * xjk;
					}
					const double djk = sqrt(djk_sq);
					// save initial end-to-end distance
					eq_dist[link_id][0] = djk;
				}
				link_id += 1;
			}
		}
	}
	void save_all_local_data(const size_t _idx)
	{
		fstream saveLink;
		saveLink.open("alllocallinkdata" + std::to_string(_idx) + ".txt", ios::app);
		saveLink << "#	j	k	lambda_loc r/n	g/n" << endl;
		for (size_t l = 0; l < Connections.size(); l++)
		{
			// computations for polymer-polymer links
			if (static_cast<size_t>(Connections[l][2]) == 0)
			{
				const size_t j = Connections[l][0];
				const size_t k = Connections[l][1];
				double djk_sq = 0.0;
				for (size_t i = 0; i < dim; i++)
				{
					const double lbox = box_size[i];
					double xjk = Nodes[j][i] - Nodes[k][i];
					xjk -= round(xjk / lbox) * lbox;
					djk_sq += xjk * xjk;
				}
				const double djk = sqrt(djk_sq);

				double seg;
				for (size_t m = 0; m < Neighbors[j].size(); m++)
				{
					if (Neighbors[j][m] == k)
					{
						seg = Segments[j][m];
						m = Neighbors[j].size();
					}
				}
				const double extension = djk / seg;
				const double lambda_local = djk / eq_dist[l][0];
				const double gn = Free_Energy[l] / seg;

				saveLink << Connections[l][0] << "		" << Connections[l][1] << "		" << lambda_local << "		" << extension << "		" << gn << endl;
			}
		}
		saveLink.close();
	}

	////////////////////////////////////////////////////////////////////////
public:
	network(const size_t _dimension)
	{
		dim = _dimension;
		box_size.resize(dim);

		Links_Node = 2 * dim;
		LinksN = 0;
		Links_per_Node = 0.0;

		NodesN = 0;
		Nodes.clear();
		init_length = 0.0;
		kappa = 0.0;
		defect_frac = 0.0;
		crystallizing = true;

		filled = false;
		filler_frac = 0.0;
		filler_bonds = Links_Node;

		rcut = 0.0;
		rcutsq = 0.0;
		rm = 0.0;
		rmsq = 0.0;

		max_strain = 1.0;
		strain_steps = 0.0;

		G = 0.0;
		eta = 0.0;
		eta_c = 0.0;

		A0_inv = 0.0;
		ratio_A1A2 = 0.0;
		sigma = 0.0;

		chi = 0.0;

		CrackCheck = false;
		RuptureCheck = false;
		CrackCentered = false;
		CrackSize = 0.0;
		RuptureEnergy = 0.0;
		RuptureCryst = 1.0;
		LinksCracked = 0;

		LinksCrackedRN = 0;

		save = false;
		// save_iteration = 0;
		save_counter = 0;

		read_initial_config = false;

		make_force_hist = false;
		Save_Segments = false;

		use_morphology_generator = false;
		no_ff = 0;
		no_pf = 0;
		MCsteps = 0;

		sigma1 = 0.0;
		sigma2 = 0.0;
	};

	void set_approx_nodes(const size_t _nodes)
	{
		NodesN = _nodes;
	}

	void set_kappa(const double _kappa)
	{
		kappa = _kappa;
		std::cout << "# Measure for perturbation of the inhomogenous lattice kappa is " << kappa << "." << endl;
	}

	void set_defects(const double _defects)
	{
		defect_frac = _defects;
	}

	void set_crystallization(const bool _crystallizing)
	{
		crystallizing = _crystallizing;
		if (crystallizing == true)
		{
			std::cout << "# The network is able to crystallize." << endl;
		}
		else
		{
			std::cout << "# The network is NOT able to crystallize." << endl;
		}
	}

	void set_eta(const double _eta)
	{
		eta = _eta;
		eta_c = eta * lc;
	}

	void set_data_saving(const bool _save, const bool _local_link_data)
	{
		if (_save == true)
		{
			save = true;
			// save_iteration = _iterations;
			if (_local_link_data == true)
			{
				local_link_data_save = true;
			}
			else
			{
				local_link_data_save = false;
			}
		}
		else
		{
			save = false;
		}
	}

	void set_cut_off(const double _cutoff)
	{
		rcut = _cutoff * r0;
		rcutsq = rcut * rcut;
		rm = rcut + 1.0;
		rmsq = rm * rm;
	}

	void set_rupture(const bool _rupture, const double _rupture_energy, const double _max_cryst_rupt)
	{
		if (_rupture == true)
		{
			RuptureCheck = true;
			RuptureEnergy = _rupture_energy;
			RuptureCryst = _max_cryst_rupt;
		}
		else
		{
			RuptureCheck = false;
		}
	}

	void set_crack(const bool _crack, const double _crack_size)
	{
		if (_crack == true && _crack_size > 0.0)
		{
			CrackCheck = true;
			CrackSize = _crack_size;
		}
		else
		{
			CrackCheck = false;
		}
	}

	void set_crack_center(const bool _crack, const double _crack_size)
	{
		if (_crack == true && _crack_size > 0.0)
		{
			CrackCheck = true;
			CrackSize = _crack_size;
			CrackCentered = true;
		}
		else
		{
			CrackCheck = false;
		}
	}

	void set_filler(const bool _filled, const double _filler_frac, const size_t _filler_bonds)
	{
		filled = _filled;
		if (_filled == true)
		{
			filler_frac = _filler_frac;
			filler_bonds = _filler_bonds;
		}
	}

	void set_morphology_generator(const bool _use_morphology_generator, const double _filler_frac, const size_t _MC_attempts)
	{
		use_morphology_generator = _use_morphology_generator;
		if (use_morphology_generator == true)
		{
			filler_frac = _filler_frac;
			MCattempts = _MC_attempts;
		}
	}

	void set_force_hist(const bool _make_force_hist, const vector<double> _force_hist_lambda)
	{
		if (_make_force_hist == true)
		{
			make_force_hist = true;
			force_hist_lambda = _force_hist_lambda;
		}
	}

	void save_Kuhn_segments(const bool _save_segments)
	{
		if (_save_segments == true)
		{
			Save_Segments = true;
		}
		else
		{
			Save_Segments = false;
		}
	}

	// this function adjusts the number of nodes, determines the corresponding length
	// of the square or cubic simulation volume, places the nodes on a sc lattice and calls the functions
	// which set up the links
	void initial_lattice(const bool _read_init_config, const std::string &_node_positions_file, const std::string &_linkage_file, const std::string &_filler_file)
	{
		// correct the number of nodes to make it fit a square or cubic lattice in 2D or 3D respectively
		const double nodes1double = pow(static_cast<double>(NodesN), static_cast<double>(1.0 / dim));

		NodesN = static_cast<size_t>(round(pow(nodes1double, static_cast<double>(dim))));
		Nodes.resize(NodesN);

		// compute total number of MC steps for the morphology generator
		// based on the number of nodes and the fixed number of MC interchange attempts per node
		if (use_morphology_generator)
		{
			MCsteps = MCattempts * NodesN;
		}

		std::cout << "# Initial number of nodes after recalculation is " << NodesN << "." << endl;

		// compute length of the simulation box
		init_length = static_cast<double>(init_dist) * nodes1double;
		for (size_t i = 0; i < dim; i++)
		{
			box_size[i] = init_length;
		}

		std::cout << "# Initial length of the simulation box is " << init_length << "." << endl;

		// compute original cross section of the sample
		A0_inv = 1.0 / pow(init_length, dim - 1.0);

		if (_read_init_config == false)
		{	
			// place nodes on a simple cubic lattice
			// in the 2D case, the 3rd dimension will be neglected when setting up the vector Nodes
			const size_t nodes1d = static_cast<size_t>(round(nodes1double));
			vector<vector<double>> init_pos;
			const double half_d_init = 0.5 * d_init;
			for (size_t h = 0; h < nodes1d; h++)
			{
				for (size_t k = 0; k < nodes1d; k++)
				{
					for (size_t l = 0; l < nodes1d; l++)
					{
						const double z = static_cast<double>(init_dist * h) + half_d_init;
						const double y = static_cast<double>(init_dist * k) + half_d_init;
						const double x = static_cast<double>(init_dist * l) + half_d_init;

						init_pos.push_back(vector<double>{x, y, z});
					}
				}
			}

			// save the initial positions in the vector Nodes, which is of size NodesNx2 or NodesNx3 for a 2D or 3D system respectively
			//savePos.open("InitPos.txt", ios::app);
			for (size_t i = 0; i < dim; i++)
			{	

				const double lbox = box_size[i];
				for (size_t j = 0; j < NodesN; j++)
				{
					// if the lattice should be inhomogenous, i.e. kappa!=0, the homogenous lattice is peturbed here
					const double perturbation = kappa * d_init * randomnumber(-1.0, 1.0);
					const double pert_coord = init_pos[j][i] + perturbation;

					// periodic boundary conditions
					init_pos[j][i] = pert_coord - round(pert_coord / lbox) * lbox;

					Nodes[j].push_back(init_pos[j][i]);
					//savePos<<init_pos[j][i]<<std::endl;
				}
			}
			//savePos.close();

			if (use_morphology_generator)
			{
				run_morphology_generator();
				nodes_neighbor_list();
				set_initial_links_morphology();
			}
			else
			{
				nodes_neighbor_list();
				set_links();
			}
		}
		else
		{	
			read_initial_config = true;
			read_initial_node_positions(_node_positions_file);
			read_initial_links(_linkage_file);
		}

		// remark: in the case of reading in the data for the initial configuration,
		// the neighbor list of links should be the same if the rest of the initial network is also identical
		links_fixed_neighbor_list_morphology();

		std::cout << "# Setting up the initial configuration is completed." << endl;
	}

	// this function starts the simulation
	void run_simulation(const double _max_strain, const double _strain_steps, const size_t _run, const size_t _cycles, const bool _stretching_only)
	{
		bool stretch_only = _stretching_only;
		if (_cycles > 1 && _stretching_only)
		{
			stretch_only = false;
		}
		max_strain = _max_strain;
		strain_steps = _strain_steps;

		size_t force_hist_index = 0;

		const int N_strains = static_cast<int>((max_strain - 1.0) / strain_steps) + 1;

		fstream results;
		results.open("results" + std::to_string(_run) + ".txt", ios::app);
		results << "# dimension = " << dim << endl;
		results << "# No of nodes = " << NodesN << endl;
		results << "# kappa = " << kappa << endl;
		results << "# eta = " << eta << endl;
		results << "# No of links = " << LinksN << endl;
		results << "# No of filler-filler links = " << no_ff << endl;
		results << "# No of polymer-filler links = " << no_pf << endl;

		if (crystallizing == true)
		{
			results << "# The network is able to crystallize." << endl;
			results << "# defect fraction = " << defect_frac << endl;
		}
		else
		{
			results << "# The network is NOT able to crystallize." << endl;
		}
		if (filled == true)
		{
			results << "# The filler faction is = " << filler_frac << " and the fillers have " << filler_bonds << " bonds." << endl;
		}
		else
		{
			results << "# There are no fillers in the network." << endl;
		}
		if (CrackCheck == true)
		{
			results << "# crack size = " << CrackSize << endl;
			if (CrackCentered == true)
			{
				results << "# The crack is inserted in the center of the network." << endl;
			}
		}
		else
		{
			results << "# No crack inserted. " << endl;
		}
		if (RuptureCheck == true)
		{
			results << "# rupture energy = " << RuptureEnergy << endl;
		}
		else
		{
			results << "# Rupture of links is not activated. " << endl;
		}

		results << "# lambda		sigma			chi			broken links		fraction of links that broke		total number of links		fraction of links with r>n			number of links with r>n that broke		filler-filler total		polymer-filler total" << endl;

		double prev_lambda = 1.0;	 // to save the previous strain for rescaling of node positions
		size_t prev_LinksN = LinksN; // total number of links at the previous strain;

		// continue preparation of the initial network
		if (use_morphology_generator == true and read_initial_config == false)
		{
			initial_FIRE_optimizer();
			set_filler_links();
		}
		// save initial configuration
		save_initial_r_of_links();
		save_data();

		for (size_t cyc = 0; cyc < _cycles; cyc++)
		{
			// save the number of segments, length, crystallinity of broken links
			LinksRuptured.open("links" + std::to_string(_run) + "cycle" + std::to_string(cyc) + ".txt", ios::app);
			LinksRuptured << "# n		r		c		c*lc		crystallinity of this link=clc/n		defect?		initial r" << endl;

			// stretch the network until the maximum strain is reached
			for (int s = 0; s < N_strains; s++)
			{
				const double lambda = 1.0 + strain_steps * static_cast<double>(s);
				cout << "# Lambda = " << lambda << endl;

				LinksRuptured << "# lambda = " << lambda << endl;

				LinksCracked = 0;
				LinksCrackedRN = 0;

				stretching(lambda, prev_lambda);

				FIRE_optimizer_morphology();

				filler_reversible_bonds();

				crystallization_check();

				FIRE_optimizer_morphology();

				filler_reversible_bonds();

				save_data();

				if (make_force_hist == true && lambda == force_hist_lambda[force_hist_index])
				{
					force_and_energy(false, true, force_hist_index);
					force_hist_index += 1;
				}

				force_and_energy(true, false, 0);
				stress();
				crystallinity_index();

				if (s == 0 && cyc == 0)
				{
					save_initial_r_of_links();
				}

				size_t LinksRN = 0;
				for (size_t m = 0; m < LinksN; m++)
				{
					if (LinkCheck[m] == true)
					{
						LinksRN += 1;
					}
				}

				const double linksChecked = static_cast<double>(LinksRN) / static_cast<double>(LinksN);

				const double link_fraction = static_cast<double>(LinksCracked) / static_cast<double>(prev_LinksN);
				results << lambda << "			" << sigma << "			" << chi << "			" << LinksCracked << "			" << link_fraction
						<< "		" << prev_LinksN << "				" << linksChecked << "				" << LinksCrackedRN
						<< "				" << no_ff << "				" << no_pf << endl;

				prev_LinksN = LinksN;

				prev_lambda = lambda;

				// exit(0);
			}

			if (stretch_only == false)
			{
				// contract the network from maximum strain to lambda = 1
				for (int s = N_strains - 2; s >= 0; --s)
				{
					const double lambda = 1.0 + strain_steps * static_cast<double>(s);
					cout << "# Lambda = " << lambda << endl;

					LinksRuptured << "# lambda = " << lambda << endl;

					LinksCracked = 0;
					LinksCrackedRN = 0;

					stretching(lambda, prev_lambda);

					FIRE_optimizer_morphology();

					filler_reversible_bonds();

					crystallization_check();

					FIRE_optimizer_morphology();

					filler_reversible_bonds();

					save_data();

					if (make_force_hist == true && lambda == force_hist_lambda[force_hist_index])
					{
						force_and_energy(false, true, force_hist_index);
						force_hist_index += 1;
					}

					force_and_energy(true, false, 0);

					stress();
					crystallinity_index();

					size_t LinksRN = 0;
					for (size_t m = 0; m < LinksN; m++)
					{
						if (LinkCheck[m] == true)
						{
							LinksRN += 1;
						}
					}

					const double linksChecked = static_cast<double>(LinksRN) / static_cast<double>(LinksN);

					const double link_fraction = static_cast<double>(LinksCracked) / static_cast<double>(prev_LinksN);

					results << lambda << "			" << sigma << "			" << chi << "			" << LinksCracked << "			" << link_fraction
							<< "		" << prev_LinksN << "				" << linksChecked << "				" << LinksCrackedRN
							<< "				" << no_ff << "				" << no_pf << endl;

					prev_LinksN = LinksN;

					prev_lambda = lambda;
				}
			}
			LinksRuptured.close();
		}
		results.close();
	}
};

#endif
