First Simulation

This tutorial provides a complete, step-by-step guide to performing your first phase separation molecular dynamics simulation using DROPPS.

Note that DROPPS is a highly flexible package that can be used in many different ways, but this tutorial focuses on the most straightforward workflow that closely follows the slab-simulation pipeline. This is the recommended starting point for new users, as it covers all the essential steps while maintaining a familiar structure.

DROPPS is intentionally designed to mimic the user interface of GROMACS, meaning that each step corresponds to a familiar GROMACS command. The following are the DROPPS equivalents of common (selected) GROMACS commands:

GROMACS	DROPPS
gmx pdb2gmx	dps pdb2dps
gmx editconf	dps editconf
gmx grompp	dps grompp
gmx mdrun	dps mdrun
gmx make_ndx	dps make_ndx
gmx contact	dps contact

This tutorial walks you through the complete pipeline for a slab simulation for a multi-component co-phase separation system (co-phase separation of α-synuclein and full-length TDP-43 ).

Overview of the Workflow

A DROPPS simulation consists of the following steps:

1. Generate single-chain structure & topology for each type of biomolecule (pdb2dps)

2. Build a multi-chain system (genmesh)

3. Prepare a simulation run input (grompp)

4. Run MD simulation on local machine or on HPC cluster (mdrun)

5. (Optional) Build index groups (make_ndx)

6. (Optional) Perform trajectory analysis (contact, density, gyrate, etc.)

Each step directly corresponds to a command in the dps command-line suite.

1. Generate Structures and Topologies (pdb2dps)

To create the initial single-chain structure and topology from a sequence or PDB file, use dps pdb2dps.

1.1 Generate Single-Chain Structure and Topology for α-Synuclein

Assume that you don’t have the all-atom structure for α-synuclein, thus the fastest way is to generate its coarse grained structure and topology is through a random walk.

dps pdb2dps
-s MDVFMKGLSKAKEGVVAAAEKTKQGVAEAAGKTKEGVLYVGSKTKEGVVHGVATVAEKTKEQVTNVGGAVVTGVTAVAQKTVEGAGSIAAATGFVKKDQLGKNEEGAPQEGILEDMPVDPDNEAYEMPSEEGYQDYEPEA
-cNTD -cCTD -r 3
-oc asyn -op asyn -on asyn -ff HPS

The command line options are described as follows:

• The -cNTD and -cCTD flags specify that additional charge is added to N- and C-terminus.

• The -oc, -op, -on flags specify the coordinate file name, topology file name, and molecule name.

• The -ff flag specify the forcefield. If not specified, you will be asked to choose interactively.

• The -r flag specifies the maximum radius for the random walk strou can adjust this value to generate a more compact or more extended initial structure.ucture. Y

The command will write the following files:

• asyn.pdb — coarse-grained coordinates

• asyn.itp — single-chain topology

Optional: Add Angle Potentials for α-Synuclein to Capture Secondary Structure Propensities

If you want to include angle potentials to capture secondary-structure propensities, you can use the dps addangle command.

If you want to include angle potentials to capture secondary-structure propensities, you can use the dps addangle command.

A recent publication:

Y. Tang, J. Wen, Z. Yang, et al., “Dissecting the Molecular Determinants of α-synuclein Phase Separation and Condensate Aging: The Pivotal Role of β-Sheet-Rich Motifs”, Advanced Science 12 (44) (2025): e11545.

showed that α-synuclein has seven β-sheet-rich domains. Their β-probabilities are listed in the following table. It is wise to add angle potentials to these domains with force constants proportional to their β-probabilities, as this may better capture the conformational ensemble of α-synuclein.

Sequence	β-probability	Angle	Force constant
4FMKGLSK10	0.83	180°	0.83
26VAEAAGK32	0.82	180°	0.82
36GVLYVGS42	0.66	180°	0.66
61EQVTNVG67	0.74	180°	0.74
75TAVAQKT81	0.77	180°	0.77
88IAAATGF94	0.66	180°	0.66
113LEDMPVD119	0.44	180°	0.44

To add these angle potentials, create a file named angle.txt with the contents provided at the end of this page. Or, use this command to create the file directly from the command line:

awk 'BEGIN{split("4 26 36 61 75 88 113",s);split("83 82 66 74 77 66 44",v);for(i=1;i<=7;i++)for(j=0;j<7;j++)printf"%d\t180\t%d\n",s[i]+j,v[i]}' > angle.txt

Then, run the following command to add angle potentials:

mv asyn.itp asyn_base.itp
dps addangle -ip asyn_base.itp -op asyn.itp -al angle.txt

This command will generate a new topology file that includes the angle potentials.

1.2 Generate Single-Chain Structure and Topology for Full-Length TDP-43

TDP-43 contains multiple structured domains whose conformation should be resumed in coarse-grained simulations. However, the all-atom structure for full-length TDP-43 is not available. To generate the coarse-grained structure and topology for full-length TDP-43, we can use the predicted structure from AlphaFold. The predicted structure can be downloaded from the AlphaFold Protein Structure Database (https://alphafold.ebi.ac.uk/entry/AF-Q9W1I0-F1).

After downloading the PDB file, you can run the following command to generate the coarse-grained structure and topology:

dps pdb2dps
-s MSEYIRVTEDENDEPIEIPSEDDGTVLLSTVTAQFPGACGLRYRNPVSQCMRGVRLVEGILHAPDAGWGNLVYVVNYPKDNKRKMDETDASSAVKVKRAVQKTSDLIVLGLPWKTTEQDLKEYFSTFGEVLMVQVKKDLKTGHSKGFGFVRFTEYETQVKVMSQRHMIDGRWCDCKLPNSKQSQDEPLRSRKVFVGRCTEDMTEDELREFFSQYGDVMDVFIPKPFRAFAFVTFADDQIAQSLCGEDLIIKGISVHISNAEPKHNSNRQLERSGRFGGNPGGFGNQGGFGNSRGGGAGLGNNQGSNMGGGMNFGAFSINPAMMAAAQAALQSSWGMMGMLASQQNQSGPSGNNQNQGNMQREPNQAFGSGNNSYSGSNSGAAIGWGSASNAGSGSGFNGGFGSSMDSKSSGWGM
-f AF-Q9W1I0-F1-model_v6.pdb
-cNTD -cCTD
-oc tdp43.pdb -op tdp43.top -on tdp43 -ff HPS

The command line options are described as follows:

• The -f flag specifies the input PDB file containing the all-atom structure.

• The -cNTD and -cCTD flags specify that additional charge is added to N- and C-terminus.

• The -oc, -op, -on flags specify the coordinate file name, topology file name, and molecule name.

• The -ff flag specify the forcefield. If not specified, you will be asked to choose interactively.

The command will write the following files:

• tdp43.pdb — coarse-grained coordinates

• tdp43.itp — single-chain topology

Optional: Add Elastic Network Bonds to Preserve Structured Domain Conformations

If you want to preserve the conformations of the structured domains in TDP-43, you can add elastic network bonds to these domains. This can be done using the dps genelastic command. For example, if you want to add elastic bonds to the NTD domain (residues 1–77), the RRM1 domain (residues 106–177) and the RRM2 domain (residues 192–259) respectively, you can create a file named elastic.txt with the content provided at the end of this page. Or, use this command to create the file directly from the command line:

seq -s' ' 1 77 > elastic.txt
seq -s' ' 106 177 >> elastic.txt
seq -s' ' 192 259 >> elastic.txt

Then, to add elastic network bonds, run the following command:

mv tdp43.itp tdp43_base.itp
dps genelastic -f tdp43.pdb -p tdp43_base.itp -o tdp43.itp -er elastic.txt

This command will generate a new topology file that includes the elastic network bonds for the specified domains. Note that the parameters for the elastic bonds can be adjusted by command line options. Please refer to the help information in “dps genelastic -h”.

2. Build a Multi-Chain System (genmesh)

To construct a system with many copies of the protein, use dps genmesh. For example, we would like to simulate a system with 125 chains of α-synuclein and 125 chains of TDP-43 on a random mesh, use the following command:

dps genmesh
-f asyn.pdb tdp43.pdb
-p asyn.itp tdp43.itp
-n 125 125
-mesh 8 8 8 -g 1 --shuffle
-oc system.pdb -op system.top

The command line options are described as follows:

• The -f flag specifies the input PDB file for each molecule specie.

• The -p flag specifies the input topology file for each molecule specie.

• The -n flag specifies the number of chains for each molecule specie.

• The -mesh flag specifies the dimensions of the spatial mesh to arrange the chains.

• The g flag specifies the minimum distance between chains to avoid overlaps.

• The –shuffle flag randomizes the placement of chains on the mesh.

• The -oc and -op flags specify the output coordinate file name and topology file name for the multi-chain system.

Note that there are more mesh points (8*8*8=512) then the number of chains (125+125=250), so the chains will be randomly placed on the mesh with some empty points.

The command will write the following files:

• system.pdb — coarse-grained coordinates for the whole system

• system.top — Topology for the whole system

3. Perform a Short NPT Simulation for Dense-Phase Formation

Before running a long production simulation, it is often helpful to perform a simulation under the NPT ensemble for a short time to allow the system form a dense phase. This is because starting from a random configuration, it may take a long time for the system to spontaneously phase separate and form a dense phase, which can be computationally expensive. By performing a short simulation under the NPT ensemble, you can allow the system to quickly phase separate and form a dense phase, which can then be used as the starting point for a longer production simulation.

To do so, you will need a molecular dynamics parameter file (npt.mdp) that specifies the simulation parameters for the NPT simulation. The content for the file and the description for each parameter is provided in the end of this page. You can adjust the parameters as needed, but for a quick test, you can use the default parameters provided in the file.

The most important parameters are:

• tcoulp = Bussi, which is the method for temperature coupling. Bussi thermostat

  is a stochastic velocity rescaling thermostat that can generate the correct canonical ensemble. It is suitable for coarse-grained simulations and is currently the only allowed method for temperature coupling in DROPPS.

• initial-temperature = 150 and production-temperature = 150 K, which is the ensemble temperature for the system.

  Use of a low temperature allow the system to efficiently collapse into a dense phase.

• pcoulp = True, which enables pressure coupling. This allows the system to adjust its

  volume and density during the simulation.

• ref-P = 1 and tau-P = 5, which is the reference pressure and time constant for

  pressure coupling. This is a common choice for simulating biological systems at ambient conditions.

Then, you can run the NPT simulation using the following two command:

dps grompp -f system.pdb -p system.top -m npt.mdp -o npt.tpr
dps mdrun -s npt.tpr -o npt

The grompp command generates a run input file (tpr) that is used for the mdrun command, which will be described extensively below.

The mdrun command will generate these files:

• npt.xtc: trajectory file

• npt.pdb: final configuration

• npt.log: log information

4. Generate a Phase-Coexistence State

After the NPT simulation, you are now in a place to prepare a phase-coexistence state. This can be done by expanding the final configuration of the NPT simulation in one direction by the following command:

dps editconf -f npt.pdb -mz 10 -o slab.pdb

The command line options are described as follows:

• The -mz flag specifies that the box will be expanded in the z direction for 10 times.

Note that the periodic boundary condition will be automatically treated so that no molecule will be break in the z direction.

5. Prepare the Production Run (Phase-Coexistence Simulation)

DROPPS uses a distributed design in which the simulation-construction and trajectory-analysis modules are decoupled from the production-run module and linked through a portable run-input file (tpr). This architecture enables simulations to be executed on external computational resources such as high-performance computing (HPC) clusters and supercomputers, with OpenMM used as backends.

For the production run, you also need a mdp file (md.mdp) that specifies the simulation parameters for the production simulation. The content for the file and the description for each parameter is provided in the end of this page. You can adjust the parameters as needed.

Note that the pressure coupling is turned off in the production simulation.

A run input file (md.tpr) can be generated using the following command:

dps grompp -m md.mdp -f system.pdb -p system.top -o md.tpr

This command will generate a “md.tpr” run input file that includes everything needed to run the production simulation.

6. Run the Production Simulation

The simulation can be performed on local machine. For better performance, it is recommended to run the simulation on an HPC cluster or supercomputer with GPU resources. The tpr file can be transferred to an HPC cluster or supercomputer for running the production simulation. A sample job submission script for a slurm-based HPC cluster is provided:

#!/bin/bash

module load DROPPS/0.3.1
dps mdrun -s md.tpr -o md >> md.out &> md.error

The mdrun command will generate these files:

• md.xtc: trajectory file

• md.pdb: final configuration

• md.log: log information

If you would like to track the progress of the simulation, or the time evolution of some critical properties such as potential energy, temperature, etc., you can adjust the settings in the mdp file section “Output control” to print these properties to the screen or the log file at a certain frequency.

Currently, these properties are allowed:

• step: current step index

• time: current time

• potentialEnergy: potential energy

• kineticEnergy: kinetic energy

• totalEnergy: total energy

• temperature: instantaneous temperature

• volume: periodic box volume

• density: system density

• progress: current progress (percent completion).

• remainingTime: estimate of the remaining clock time until completion.

• speed: estimate of the simulation speed in ns/day

• elapsedTime: elapsed time of the simulation in seconds.

7. Generate Index Groups (Optional, make_ndx)

Some analyses require custom index groups.

Use:

dps make_ndx -f md.tpr -o md.ndx

The interactive interface supports:

• splitch — automatically split groups by chain (specialized for large multichain HPS systems)

• res, abbr, chain, index, etc.

8. Perform Data Analysis (Optional)

DROPPS provides a suite of built-in analysis tools, such as:

• dps density

• dps cmap

• dps idist / dps odist

• dps gyrate

• dps contact

Example: calculate a contact number map between protein chains:

dps contact -s md.tpr -f md.xtc -n md.ndx -o contact.xvg

Example: compute density profile along the z-axis:

dps density -s md.tpr -f md.xtc -n md.ndx --axis z -o density.xvg

Troubleshooting

Common issues (summarized from user guidance on pp. 31–33):

• Command not found Ensure dps is added to your $PATH.

• Missing topology files Check that .itp and .top files are correctly generated in earlier steps.

• MDAnalysis cannot load trajectory Use dps check to diagnose tpr / xtc issues.

• Incorrect box size after genmesh Use dps editconf to adjust simulation boundaries.

Summary

You have now completed your first full DROPPS simulation, following the exact pipeline recommended by the official documentation slides. The design philosophy, commands, and workflow mirror GROMACS but are optimized for multi-chain HPS coarse-grained simulations.

Next steps:

• Learn advanced mdp parameters

• Enable elastic restraints (dps genelastic)

• Add angle potentials (dps addangle)

• Write custom analysis scripts using the trajectory API

Related files

Angle potentials for α-synuclein (angle.txt):

4    180     83
5    180     83
6    180     83
7    180     83
8    180     83
9    180     83
10   180     83
26   180     82
27   180     82
28   180     82
29   180     82
30   180     82
31   180     82
32   180     82
36   180     66
37   180     66
38   180     66
39   180     66
40   180     66
41   180     66
42   180     66
61   180     74
62   180     74
63   180     74
64   180     74
65   180     74
66   180     74
67   180     74
75   180     77
76   180     77
77   180     77
78   180     77
79   180     77
80   180     77
81   180     77
88   180     66
89   180     66
90   180     66
91   180     66
92   180     66
93   180     66
94   180     66
113  180     44
114  180     44
115  180     44
116  180     44
117  180     44
118  180     44
119  180     44

Elastic network profile for TDP-43 (elastic.txt):

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77

106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177

192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259

Molecular Dynamics Parameter File for NPT Simulation (npt.mdp)

# Integration
integrator  = Langevin       # Allowed values: Langevin, steep
friction    = 1.0            # Friction in picosecond ^ -1
dt          = 0.01           # Integration time step in ps
nsteps      = 20000000       # Number of steps to integrate
minimize    = True           # Perform constant volume minimization before this simulation
max_step    = 100000         # Maximum number of steps to perform energy minimization
seed        = 1215           # Seed for starting the simulation

comm-mode   = Linear         # Allowed values: Linear (Remove center of mass translational velocity), none.
nstcomm     = 100            # interval for center of mass motion removal

bondtype    = bond           # Treat bond as: bond or constraint

# Energy minimization
# energytol       = 10.0      # The minimization is converged when the maximum energy is smaller than this value
forcetol        = 100         # The minimization is converged when the maximum force is smaller than this value

# Output control
nst-xout     = 10000           # Number of steps that elapse between writing coordinates to the output trajectory file
nst-cp       = 0          # Number of steps that elapse between writing checkpoint

nst-screenlog = 10000  # Number of steps that elapse between writing properties to screen, 0 for write nothing.
nst-filelog = 10000   # Number of steps that elapse between writing properties to logfile, 0 for write nothing.
screenlog-grps = step,volume,elapsedTime,remainingTime,speed,progress   # Write these properties to log file.
filelog-grps = step,volume,potentialEnergy,kineticEnergy      # Write these properties to log file.

# nstvout     = 10000          # Number of steps that elapse between writing velocity to the output trajectory file

# Neighbor searching
nstlist     = 10       # Frequency to update the neighbor list

#cutoff-scheme-lj    = static  # dynamic (lj-multi   * sigma) or static (cutoff-lj)
#cutoff-scheme-coul  = static  # dynamic (coul-multi * sigma) or static (cutoff-coul)
#cutoff-lj-multi = 3.0         # Appliable only when cutoff-scheme-lj = dynamic
#cutoff-coul-multi = 3.0       # Appliable only when cutoff-scheme-coul = dynamic
cutoff-lj   = 1.5      # Default cutoff value for LJ interaction,      only when cutoff-scheme-lj = static
cutoff-coul = 1.5      # Default cutoff value for Coulomb interaction, only when cutoff-scheme-coul = static
buffer      = 0.5      # Buffer for neighbor searching

# Force calculations
coulombtype = yukawa   # Method to calculate Coulomb interaction. Allowed values: yukawa,no
vdwtype     = pLJ      # Method to calculate vdw interactions. Allowed values: pLJ,none pLJ is perturbed LJ
salt_conc   = 0.1      # Salt concentration in mol/L

# Temperature coupling
tcoulp                  = Bussi    # Method to perform temperature coupling. Allowed values: Bussi,no. Currently only Bussi is allowed.
production-temperature  = 150      # The ensemble temperature for the system.
gen-vel                 = True     # Randomly generate velocity according to init-temperature
initial-temperature     = 150      # The ensemble temperature for the system.
warming-speed           = 1        # Speed of temperature climbing if production temperature > initial temperature, in K/ns

# Pressure coupling
pcoulp   = True     # Whether to perform pressure coupling
ref-P    = 1        # The reference pressure for coupling, in bar.
tau-P    = 5        # Frequency at which Monte Carlo pressure changes should be attempted (in time steps)

Molecular Dynamics Parameter File for Production Simulation (md.mdp)

# Template parameter file for MD simulation

# Integration
integrator  = Langevin       # Allowed values: Langevin, steep
friction    = 1.0            # Friction in picosecond ^ -1
dt          = 0.02           # Integration time step in ps
nsteps      = 250000000      # Number of steps to integrate
minimize    = True           # Perform constant volume minimization before this simulation
max_step    = 100000         # Maximum number of steps to perform energy minimization
seed        = 1215           # Seed for starting the simulation

comm-mode   = Linear         # Allowed values: Linear (Remove center of mass translational velocity), none.
nstcomm     = 100            # interval for center of mass motion removal

# Energy minimization
# energytol       = 10.0      # The minimization is converged when the maximum energy is smaller than this value
forcetol        = 100         # The minimization is converged when the maximum force is smaller than this value

# Output control
nst-xout     = 10000          # Number of steps that elapse between writing coordinates to the output trajectory file
nst_cp       = 10000          # Number of steps that elapse between writing checkpoint

nst-screenlog = 1000000  # Number of steps that elapse between writing properties to screen, 0 for write nothing.
nst-filelog = 1000000   # Number of steps that elapse between writing properties to logfile, 0 for write nothing.
screenlog-grps = step,elapsedTime,remainingTime,speed,progress   # Write these properties to log file.
filelog-grps = step,potentialEnergy,kineticEnergy      # Write these properties to log file.

# nstvout     = 10000          # Number of steps that elapse between writing velocity to the output trajectory file

# Neighbor searching
#nstlist     = 10       # Frequency to update the neighbor list

#cutoff-scheme-lj    = static  # dynamic (lj-multi   * sigma) or static (cutoff-lj)
#cutoff-scheme-coul  = static  # dynamic (coul-multi * sigma) or static (cutoff-coul)
#cutoff-lj-multi = 3.0         # Appliable only when cutoff-scheme-lj = dynamic
#cutoff-coul-multi = 3.0       # Appliable only when cutoff-scheme-coul = dynamic
cutoff-lj   = 2.0      # Default cutoff value for LJ interaction,      only when cutoff-scheme-lj = static
cutoff-coul = 2.0      # Default cutoff value for Coulomb interaction, only when cutoff-scheme-coul = static
#buffer      = 0.5      # Buffer for neighbor searching

# Force calculations
coulombtype = yukawa   # Method to calculate Coulomb interaction. Allowed values: yukawa,no
vdwtype     = pLJ      # Method to calculate vdw interactions. Allowed values: pLJ,none pLJ is perturbed LJ
salt_conc   = 0.1      # Salt concentration in mol/L

# Temperature coupling
tcoulp                  = Bussi    # Method to perform temperature coupling. Allowed values: Bussi,no. Currently only Bussi is allowed.
production-temperature  = 310      # The ensemble temperature for the system.
gen-vel                 = True     # Randomly generate velocity according to init-temperature
initial-temperature     = 150      # The ensemble temperature for the system.
warming-speed           = 1        # Speed of temperature climbing if production temperature > initial temperature, in K/ns

# Pressure coupling
pcoulp   = False     # Whether to perform pressure coupling
#ref-P    = 1        # The reference pressure for coupling, in bar.
#tau-P    = 20       # Frequency at which Monte Carlo pressure changes should be attempted (in time steps)