README file for data and code associated with the "Surfactancy in a tadpole model of proteins" manuscript.

//////////////////////////////////////////////////////////////////////////
1: Running the Polywag code

The simulation code (named Polywag) is found in the polywag_DG folder.
It is comprised of many files, but only the polywag_input.c file - which sets the simulation parameters - is ever edited.
The rest specify the functions used in the simulation, and were unchanged throught the collection of data in this work.


1.1: Setting parameters in polywag_input.c

The first function in this file (get_directoty_and_suffix) sets where the files are to be read/written, and the part of the filenames consistent across the simulation.
You might want to modify the directory, set to the current location (./) by default.

set_system_parameters is where the important simulation parameters are set.
Most of these were fixed to their value in the file throughout the work. Those will not be discussed here.
The parameters that were modified were:
	- params->L[*]: this sets the box dimensions. Was {27,4,4} for most of this work.
	- params->Ncolloid[0]: sets the number of colloid particles.
	- params->Ntadpole[0]: sets the number of tadpole molecules.
	- params->Ntail[0]: sets the number of beads in the tadpole tail.
	- params->Npoly[0]: sets the number of polymer molecules.
	- params->Nmono[0]: sets the number of beads in each polymer chain.

We also had 2 pairs of potentials. For most of the work we set the following in the loop under "Parameters associated with interaction potentials":
	params->EV_fn[i][j] = &EV_pseudoHS;
        strcpy(potl_params[i][j].EV_type,"Pseudo HS");
        potl_params[i][j].EV_energy = 1.5*params->kT;
        potl_params[i][j].EV_length = 0.98 * params->diameter_ij[i][j];

	params->bond_fn[i][j] = &bond_well_quad;
        strcpy(potl_params[i][j].bond_type,"well_quad");
        potl_params[i][j].spring_const = 8.819 / (params->diameter_ij[i][j] * params->diameter_ij[i][j]);
        potl_params[i][j].bond_length = params->diameter_ij[i][j] + params->diameter_ij[i][j];

However, for the interface tension simulations we instead used:
	params->EV_fn[i][j] = &EV_WCA;
        strcpy(potl_params[i][j].EV_type,"WCA");
        potl_params[i][j].EV_energy = 1.0*params->kT;
        potl_params[i][j].EV_length = 1.0*params->diameter_ij[i][j];

	params->bond_fn[i][j] = &bond_FENE;
        strcpy(potl_params[i][j].bond_type,"FENE");
        if(i==2 && j==2){potl_params[i][j].spring_const = 7.8221296;}
        else{potl_params[i][j].spring_const = 2.5577892;}
        potl_params[i][j].bond_length = params->diameter_ij[i][j] + 2.0*params->diameter_ij[2][2];

The following output functions were set to true: system_params, species, unwrapped_config and cont_config_out.
Where we were continuing simulations across multiple jobs or starting from an pre-initialised state cont_config_in was also true.
The rest were all set to false in this work.


1.2: Running from a bash shell

The makefile contains the commands for running the Polywag code from a bash shell.
It can be compiled with the gcc compiler with the command "make polywag" or with an intel compiler using "make polywagicc".
Note the intel compiler is typically faster, but needs the appropriate module to be loaded first.

Once compiled, "make run" with start the simulation.

As seen in the run command definition, several variables are passed in as arguments. These work as follows:
	- Batch: sets the second to last number in file names. Is used to keep track of which batch jobs simulations were part of.
		(Note this is replaced by the seed if Batch == 0)
	- Index: sets the last number in file names. Is used to track simulations within a batch job.
	- Continue: specifies which batch to continue from. The data will be appended to the appropriate files if they exist.
		(Note this won't have any effect if Continue == 0 or cont_config_in == false)
	- Param: is passed into set_system_params and can be used to vary a parameter(s) externally.


//////////////////////////////////////////////////////////////////////////
2: Interpreting data files

2.1: File names

We formatted file names like so:
	*_Nc(Ncolloid[0])_Nh(Ntadpole[0])_Nt(Ntail[0])_Np(Npoly[0])_Nm(Nmono[0])_V(L[0]*L[1]*L[2])_(Batch)_(Index).txt

This quickly shows the particle numbers and box volume, as well as unique indentifiers at the end of the files.
All files with the same parameters in (...) above but different * are associated with the same simulation.


2.2: 'unwrapped' files

These are the raw data used. 
They list snapshots of particle centres with the format:
	time1
	x1	y1	z1
	x2	y2	z2
	...
	xNtot	yNtot	zNtot	
	time2
	x1	y1	z1
	x2	y2	z2
	...
	xNtot	yNtot	zNtot
etc.

Note the particles are ordered so colloids come first, then tadpoles (who occupy Ntail[0]+1 lines), then polymers (who occupy Nmono[0] lines).


2.3: 'system' files

These list system parameters in a format that is easy for humans to read.


2.4: 'species' files

These list the species numbers associated with each particle so you know which parameters they used from the 'system' file.
It also lists their molecule numbers and whether they are bonded to the particles above and below them in the list.


2.5: 'colRanges' and 'polyRanges'

These files specify ranges of positions along the long box axis where sufficiently far from the interfaces to be safely classified as bulk phases.
Each row is specific to the simulation with corresponding index within the folder.

These were used when calculating densities for the end points of tie lines in Fig. 6.
They were also used to correct for interface drift and re-centre simulations so we could average over simulations in Fig. 4.


2.6: Non-raw files

The remaining files included were produced during analysis of the 'unwrapped' data.
These are:

'g2' files list the radial distribution functions between all species, at each time.
These were created using the RDF_analysis_v2.c code.

'tension' files list the interface tension measured at each time.
These are specific to the WCA + FENE sims, and were created using the slab_analysis.c code.


//////////////////////////////////////////////////////////////////////////
3: Analysis codes

3.1: slab_analysis.c

After setting parameters, this code loops over snapshots in the appropriate 'unwrapped' files.
From there it bins the system into slabs perpendicular to the box axes.

In this way it obtains densities and stress tensors a function of position.
These are ouputted in 'x/y/zslabs' files (each with slabs perp to the x/y/z directions.
These files are not needed for our figures, and hence not included in the data provided.

For the WCA_FENE_1266 data, it also produced the 'tension' files which log the measured interface tension in each snapshot.

This code can be run using the 'slab' command defined in the makefile.


3.2: RDF_analysis_v2.c

The start of this code is similar to the slab_analysis code.
Its job is to calculate the radial distribution functions, listed in the 'g2' files, which specify the time of the snapshot above the arrays of data.
These files are used in checking systems have equilibrated.


3.3: Figs_for_tadpole_ms.nb

This Mathematica notebook contains the code used in producing the figures.
The cell in section 0 will need running to define some functions and to set the directory you have put the data in, but otherwise the sections are self-standing.

Section 1 has code producing all 6 figures, and should work without modification provided 'baseDirectory' is correct and the sub-folders for the data are unchanged.

Section 2 has code used in analysing the data, but wasn't used directly to produce figures.