This C program uses the Metropolis Monte Carlo method and dissipative particle dynamics (DPD) to investigate conformational averages and energy profiles of a polymer in a solvent.

To get latest version:

git clone https://github.com/tdunn19/dpd

To compile and run:

cd dpd/src
make opt3
./dpd.run

Edit src/dpd.inp to change parameters. The following is a description of all the input parameters as well as typical values in brackets:

n_mon - the number of monomers
density_s - the solvent density (3)
calc_list - 1 to use cell list method, 0 to use brute force method of calculation

length_x - system box length in the x (10)
length_y - in the y (10)
length_z - in the z (long enough to fit the polymer)

density_w - density of wall particles (3)
n_layers - number of wall layers (2)

pore_radius - half side length of the rectangular pore (1.0)

r_c - cutoff radius of particle interaction (1.0)
dr_max_dpd - max displacement of a dpd monte carlo move (0.2)
dr_max_mon - max displacement of a monomer monte carlo move (0.1)

a_mm - monomer-monomer interaction strength (25.0)
a_ms - monomer-solvent interaction strength (15.0)
a_ss - solvent-solvent interaction strength (25.0)
a_sw - solvent-wall interaction strength (9.01)

pol_init_z - initial z coordinate of the first monomer (<length_z/2)
pol_init_bl - initial bond length of the polymer (0.5-1.0)

n_steps - number of monte carlo moves to attempt (1000000)
mc_ratio - percentage of mc moves that will choose a monomer over a solvent particle (0.5)
temp - system temperature (1.0)
freq_sample - sample the system every x steps (100)
freq_monitor - monitor the system every x steps (100)

iseed - random seed (change this before every run)

The run.sh and submit.sh scripts were written specifically for ACEnet clusters.

Version 4.0 (June 17, 2013)

• the nanopore has been added:
  • new input parameter pore_radius
  • new parameter r_pore: distance between pore particles
  • new parameter n_pore: number of pore particles
  • new parameter pore_volume
  • new parameter n_pore_1d: 3d vector that stores number of pore particles in each dimension
  • new parameters pore_max and pore_min: vectors that define the bounds of the pore
  • system length is now adjusted to accomodate the pore while still keeping the wall density and spacing consistent
• changed function init_param: now calculates all system parameters and allocates memory for the particle arrays
• new function init_pore: places the pore wall particles (no layers)
• changed function init_wall: initializes wall on a skewed grid according to pore parameters, and adds layers to the nanopore
• new function check_pore: takes in particle position vector r and checks if it is in the pore
  • adjusts new parameter pore_overlap if the particle is in the pore wall
  • accept_move will now check for this overlap and immediately reject the move if 1
  • returns 1 if the particle is in the pore, 0 otherwise
• changed functions calc_pressure, calc_energy_dpd, calc_energy_mon and calc_energy_brute to take into account pore particles
• renamed n_dpd to n_solvent: n_dpd is now a sum of n_solvent, n_wall and n_pore
• renamed random_move_mon and random_move_dpd to move_mon and move_solvent
• some calculations were incorrectly casted as integers: e.g. (int) a / b should be (int) (a / b)

Version 3.0 (June 7, 2013)

• the wall has been added
  • new input parameter density_w: set density of wall particles
  • new input parameter n_layers: set number of wall layers
  • new parameter n_wall: number of wall particles
  • new parameter a_sw: interaction strength between solvent and wall particles
  • new parameter r_wall: distance between wall particles
• new function init_wall: calculates n_wall, allocates memory for both the wall and the solvent particles in the array part_dpd and initializes the wall particle positions
• renamed function setup_coords to init_part: initializes monomers and solvent particles
• new function check_wall: takes in a particle position vector r and checks if it is inside the wall
  • adjusts new parameter wall_overlap accordingly
  • the random move functions will now check for this overlap and immediately reject the move if necessary
• new function calc_cm: calculates center of mass of the polymer
• new function calc_rg: calculates radius of gyration of the polymer
• new input parameter mc_ratio: allows the user to adjust the frequency that monomers or solvent particles are chosen
  • values in the range 0 to 1
  • 0 = solvent particles moved only
  • 1 = monomers moved only
• new input parameter pol_init_z: sets the initial position of the first monomer of the chain
• new input parameter pol_init_bl: sets the initial bond length of the polymer
• split up sys.length into length.x, length.y and length.z to allow variable box length in each dimension
  • split up n_cell and r_cell as well to accomodate this change
• system volume calculation now subtracts the volume occupied by the wall
• modified the sample function to find absolute monomer positions (in the absence of periodic boundary conditions)
  • calc_re, calc_bond_length and calc_cm should now work properly if a periodic boundary is encountered
• renamed function periodic_bc to periodic_bc_dr: takes a displacement vector dr and adjusts according to periodic boundary conditions
• new function periodic_bc_r: takes a position vector r and adjusts according to periodic boundary conditions
• calc_pressure now accounts for wall-solvent and wall-monomer ineractions
• changed functions in init.c to store initial positions in part.r as well as part.ro
• removed a redundant call of new_list in main.c, as it is called in init.c already

Version 2.3 (May 31, 2013)

• bug fixes
• when old neighbour cells were considered, the head of chain was set to -1 to indicate that cell has already been counted. This caused calc_energy_dpd and calc_energy_mon to also not consider those cells, which was not intended
  • made a new 3D array to deal with this: sys.hoc_copy
  • can freely alter sys.hoc_copy under random_move functions without affecting calc_energy functions
  • in the future, may use memcpy function for efficiency, but getting segfaults with it now
• if a move is accepted, particles will now overwrite the old energies with the new energies
  • that is, part.Eo = part.E
  • this ensures that if a move gets rejected and all the particle energies are reset, only those affected by the move will be reset to part.Eo
• typo on lines mc.c:124-126
  • should be part_dpd, not part_mon
• rewrote periodic_bc to use pointers
• typo in calc_bond_length: i < sys.n_mon-2 should be i < sys.n_mon-1
• new input parameters:
  • dr_max_dpd and dr_max_mon: set maximum displacement of a MC move
  • n_attempt_dpd and n_attempt_mon: counters for number of MC moves
  • n_accept_dpd and n_attempt_mon: counters for number of accepted MC moves
• print_stats now prints selection and acceptance statistics of the monte carlo moves

Version 2.2 (May 27, 2013)

• added monitor structure:
  • new variables: mon.energy, mon.re2, mon.rex, mon.rey, mon.rez, mon.bond_length, monitor_step
  • new function monitor_mem: to allocate memory for monitored quantities
  • new input parameter freq_monitor: note that this must be equal to freq_sample (may be changed in the future)
• added some new sample quantities: energy, re2, re2x, re2y, re2z
  • new functions: calc_pressure, calc_re, calc_bond_length
  • calc_pressure was modified to account for polymer
• new function periodic_bc: takes in dr vector and adjusts it based on periodic boundary conditions
• new module energy.c: functions from calc.c and cell.c which related to energy were moved here
• new function calc_energy_brute: brute force calculation of all monomer and dpd particle energies (deleted old calc_energy)
• new function check_bond(i): checks for breakage in the bonds on either side of monomer i and adjust sys.bond_break accordingly
  • removed this same functionality from energy_fene, as it is now redundant
• fixed a bug: conservative soft pair potentials were missing a factor of 1/2
• renamed output to print_stats, a more descriptive name
• adjusted write_mon to include new monitored quantities
• fixed a bug: under monte_carlo, if a move was rejected, only one type of particle energy was reset (i.e. just dpd or just mon)
• still to do: change calc_re to take into account periodic boundary conditions (might be tricky)

Version 2.1 (May 25, 2013)

• some bug fixes and new functionality
• rewrote the criteria for calculating new energies after a particle moves
  • check for new cell
  • find dix, diy, diz (difference between cell i and j)
  • using these, calculate energies for old neighboring cells
  • old method: 555 = 125 cells, every time
  • new method: 333 = 27 cells, if particle stays in the same cell
  • = 27 + 9 cells, if one cell dimension changes (e.g. ix != ixo)
  • = 27 + 9 + 6 cells, if two cell dimensions change (e.g. ix, iy != ixo, iyo)
  • = 27 + 9 + 6 + 4 cells, if three cell dimensions change (e.g. ix, iy, iz != ixo, iyo, izo)
• new function accept_move: returns 1 (true) if MC move is accepted, or 0 (false) if it is rejected
• new function check_cell: returns 1 (true) if a particle has entered a new cell, 0 (false) otherwise
• periodic boundary conditions for vector dr were incorrectly incrementing/decrementing by sys.length/2 instead of sys.length
• monomers were not being initialized in the middle of the x-y plane

Version 2.0 (May 24, 2013)

• the polymer has been added to the system
• note: brute force method currently does not work with polymer
• new variables: n_mon and a_ij broken up into a_mm, a_ms and a_ss
• calc_energy_list has been split into calc_energy_dpd and calc_energy_mon
• random_move has been split into random_move_dpd and random_move_mon
• energy_ij has been renamed to energy_c for conservative
• new function energy_fene: takes in two bonded monomers, outputs FENE bond potential
  • r_max: max distance between bonded monomers before the bond breaks
  • r_eq: equilibrium bond length
  • r_0: defined as r_max - r_eq
  • bond_break: system parameter assigned the value of 1 (true) if a monomer-monomer FENE bond breaks (r_ij > r_max)

Version 1.1 (May 23, 2013)

• some efficiency changes:
• a new cell list will only be generated when a particles moves into a new cell
• new energies will only be calculated for particles in neighboring (2 neighbors deep) cells of the randomly chosen particle

Version 1.0 (May 20, 2013)

• employed cell list method of calculation

• add the Q coordinate and binning
• add functionality for a circular pore
• add functionality to completely remove the wall and/or polymer using just input parameters
• instead of using a global 3d array hoc_copy, use the memcpy function