04 - Amber simulations¶
Prerequisite(s): 03-tleap
The Amber simulation protocol transforms the solvated, parameterized system into a production-ready molecular dynamics trajectory through systematic energy minimization, thermal equilibration, and density relaxation. This multi-stage approach ensures proper system conditioning while avoiding common simulation artifacts such as structural instabilities, unrealistic conformational changes, or thermodynamic inconsistencies that can compromise the physical validity of the results.
Simulation Strategy Overview¶
The complete simulation protocol consists of three distinct phases, each addressing specific aspects of system preparation and equilibration. The minimization phase removes steric clashes and optimizes unfavorable geometries introduced during system construction. The relaxation phase gradually brings the system to target temperature and pressure conditions while allowing structural adaptation to the simulation environment. The production phase collects equilibrium trajectory data under stable thermodynamic conditions for subsequent analysis. Each phase employs carefully designed restraint strategies, thermodynamic protocols, and sampling parameters to ensure smooth transitions between stages while maintaining structural integrity of critical protein features. The progressive nature of this approach prevents sudden perturbations that could drive the system into non-physical configurations or cause simulation instabilities.
Minimization¶
The energy minimization phase systematically relaxes the initial structure through a series of constrained optimization procedures. Each minimization stage employs progressively reduced restraints, allowing the system to adapt gradually while preserving essential structural features. This approach prevents large-scale structural rearrangements that could compromise the integrity of the roGFP2 fold or disrupt critical interactions in the chromophore environment.
01_min¶
The initial minimization stage applies strong positional restraints to all non-hydrogen atoms, allowing optimization of hydrogen positions and relief of immediate steric clashes without significant structural rearrangement. This conservative approach protects the experimentally determined protein structure while addressing force field inconsistencies and geometric irregularities.
YAML
- stage_name: 01_min
topo_path: "$SLURM_SUBMIT_DIR/inputs/mol.prmtop"
prev_restart_path: "$SLURM_SUBMIT_DIR/inputs/mol.inpcrd"
ref_coord_path: "$SLURM_SUBMIT_DIR/inputs/mol.inpcrd"
input_kwargs:
imin: 1
ntx: 1
irest: 0
ntmin: 1
maxcyc: 5000
ncyc: 1000
ntr: 1
restraint_wt: 5.0
restraintmask: "!(@H=)"
ntb: 1
ntf: 1
ntc: 1
cut: 10.0
ntxo: 2
ntwr: 200
ntpr: 1
ntwx: 200
ioutfm: 1
iwrap: 1
The restraint weight of 5.0 kcal/mol/Ų provides sufficient constraint to maintain the overall protein architecture while allowing local relaxation.
The restraint mask !(@H=)
targets all atoms except hydrogens, ensuring that hydrogen positions can optimize freely to establish proper hydrogen bonding networks.
The 5000-step limit with 1000 steepest descent steps followed by conjugate gradient optimization provides adequate convergence for this constrained system.
02_min¶
The second minimization stage relaxes restraints to exclude water molecules and ions, allowing the solvent to optimize around the restrained protein structure. This selective approach enables proper solvation shell formation while maintaining protein integrity and preventing unrealistic structural distortions.
YAML
- stage_name: 02_min
topo_path: "$SLURM_SUBMIT_DIR/inputs/mol.prmtop"
ref_coord_path: "$SLURM_SUBMIT_DIR/inputs/mol.inpcrd"
input_kwargs:
imin: 1
ntx: 1
irest: 0
ntmin: 1
maxcyc: 5000
ncyc: 1000
ntr: 1
restraint_wt: 5.0
restraintmask: "!(:WAT) & !(@H=) & !(:Na+,Cl-)"
ntb: 1
ntf: 1
ntc: 1
cut: 10.0
ntxo: 2
ntwr: 200
ntpr: 1
ntwx: 200
ioutfm: 1
iwrap: 1
The modified restraint mask !(:WAT) & !(@H=) & !(:Na+,Cl-)
excludes water molecules, hydrogens, and ions from positional restraints while maintaining constraints on protein heavy atoms.
This approach allows the solvent and electrolyte to equilibrate around the protein while preserving the structural integrity of the folded state.
03_min¶
The third minimization stage transitions to backbone-only restraints, allowing side chain relaxation while maintaining secondary structure integrity. This approach recognizes that side chain conformations may require significant adjustment to accommodate the simulation force field while preserving the essential beta-barrel architecture of the roGFP2 fold.
YAML
- stage_name: 03_min
topo_path: "$SLURM_SUBMIT_DIR/inputs/mol.prmtop"
ref_coord_path: "$SLURM_SUBMIT_DIR/inputs/mol.inpcrd"
input_kwargs:
imin: 1
ntx: 1
irest: 0
ntmin: 1
maxcyc: 5000
ncyc: 1000
ntr: 1
restraint_wt: 2.0
restraintmask: "!(:WAT) & (@C,CA,N,O,O5',P,O3',C3',C4',C5')"
ntb: 1
ntf: 1
ntc: 1
cut: 10.0
ntxo: 2
ntwr: 200
ntpr: 1
ntwx: 200
ioutfm: 1
iwrap: 1
The reduced restraint weight of 2.0 kcal/mol/Ų and backbone-specific mask (@C,CA,N,O,O5',P,O3',C3',C4',C5')
allow greater structural flexibility while maintaining the protein fold.
This restraint pattern targets the polypeptide backbone and nucleic acid structural elements, providing a foundation for side chain optimization.
04_min¶
The final minimization stage further reduces restraint strength while maintaining backbone constraints, allowing near-complete structural relaxation while preserving fold integrity. This stage provides the optimized starting structure for subsequent thermal equilibration.
YAML
- stage_name: 04_min
topo_path: "$SLURM_SUBMIT_DIR/inputs/mol.prmtop"
ref_coord_path: "$SLURM_SUBMIT_DIR/inputs/mol.inpcrd"
input_kwargs:
imin: 1
ntx: 1
irest: 0
ntmin: 1
maxcyc: 5000
ncyc: 1000
ntr: 1
restraint_wt: 1.0
restraintmask: "!(:WAT) & (@C,CA,N,O,O5',P,O3',C3',C4',C5')"
ntb: 1
ntf: 1
ntc: 1
cut: 10.0
ntxo: 2
ntwr: 200
ntpr: 1
ntwx: 200
ioutfm: 1
iwrap: 1
The minimal restraint weight of 1.0 kcal/mol/Ų provides gentle guidance while allowing substantial conformational flexibility. This final optimization ensures that the system reaches a stable energy minimum before beginning molecular dynamics simulation.
Relaxation¶
The equilibration phase gradually brings the minimized system to target simulation conditions through controlled heating and pressure equilibration. This process requires careful management of temperature and pressure coupling to avoid system instabilities while ensuring proper thermodynamic equilibration.
05_relax_nvt_r¶
The heating phase employs constant volume (NVT) conditions to gradually raise the system temperature from 100 K to 300 K over 20 ps of simulation time. This controlled heating prevents thermal shock while allowing kinetic energy to distribute throughout the system. Backbone restraints maintain structural integrity during the heating process.
YAML
- stage_name: 05_relax_nvt_r
topo_path: "$SLURM_SUBMIT_DIR/inputs/mol.prmtop"
prev_restart_path: "$SLURM_SUBMIT_DIR/../../03-min/outputs/04_min.rst"
ref_coord_path: "$SLURM_SUBMIT_DIR/inputs/mol.inpcrd"
input_kwargs:
irest: 0
ntx: 1
ig: -1
dt: 0.002
nstlim: 10000
nscm: 500
ntr: 1
restraint_wt: 1.0
restraintmask: "!(:WAT) & (@C,CA,N,O,O5',P,O3',C3',C4',C5')"
ntb: 1
ntf: 2
ntc: 2
cut: 10.0
ntt: 3
tempi: 100.0
temp0: 300.0
gamma_ln: 5.0
ntp: 0
ntxo: 2
ntwr: 500
ntpr: 500
ntwx: 500
ioutfm: 1
iwrap: 1
The Langevin thermostat (ntt: 3
) with a collision frequency of 5.0 ps⁻¹ provides efficient temperature control while introducing appropriate stochastic forces.
The initial temperature of 100 K and target temperature of 300 K create a controlled heating ramp over the 10,000-step simulation.
The backbone restraint weight of 1.0 kcal/mol/Ų maintains structural stability during thermal activation.
06_relax_npt_r¶
The pressure equilibration phase transitions to NPT conditions, allowing the system volume to adjust to the target pressure while maintaining temperature control. Reduced restraints permit greater structural flexibility while ensuring continued fold stability during density equilibration.
YAML
- stage_name: 06_relax_npt_r
topo_path: "$SLURM_SUBMIT_DIR/inputs/mol.prmtop"
ref_coord_path: "$SLURM_SUBMIT_DIR/inputs/mol.inpcrd"
input_kwargs:
irest: 1
ntx: 5
ig: -1
dt: 0.002
nstlim: 500000
nscm: 500
ntr: 1
restraint_wt: 0.5
restraintmask: "!(:WAT) & (@C,CA,N,O,O5',P,O3',C3',C4',C5')"
ntb: 2
ntf: 2
ntc: 2
cut: 10.0
ntt: 3
temp0: 300.0
gamma_ln: 5.0
ntp: 1
barostat: 2
pres0: 1.01325
mcbarint: 100
comp: 44.6
taup: 1.0
ntxo: 2
ntwr: 5000
ntpr: 500
ntwx: 5000
ioutfm: 1
iwrap: 1
The extended 1 ns simulation (nstlim: 500000
) allows sufficient time for volume equilibration under the Monte Carlo barostat (barostat: 2
).
The reduced restraint weight of 0.5 kcal/mol/Ų permits greater conformational flexibility while maintaining essential structural features.
The compressibility of 44.6 × 10⁻⁶ bar⁻¹ approximates the experimental value for liquid water.
07_relax_npt¶
The final equilibration stage removes all positional restraints, allowing the protein to adopt its preferred conformation under the simulation force field and environmental conditions. This unrestrained equilibration ensures that the system reaches a proper thermodynamic equilibrium before production data collection.
YAML
- stage_name: 07_relax_npt
topo_path: "$SLURM_SUBMIT_DIR/inputs/mol.prmtop"
ref_coord_path: "$SLURM_SUBMIT_DIR/inputs/mol.inpcrd"
input_kwargs:
irest: 1
ntx: 5
ig: -1
dt: 0.002
nstlim: 500000
nscm: 500
ntr: 0
ntb: 2
ntf: 2
ntc: 2
cut: 10.0
ntt: 3
temp0: 300.0
gamma_ln: 5.0
ntp: 1
barostat: 2
pres0: 1.01325
mcbarint: 100
comp: 44.6
taup: 1.0
ntxo: 2
ntwr: 5000
ntpr: 500
ntwx: 5000
ioutfm: 1
iwrap: 1
The absence of restraints (ntr: 0
) allows complete conformational freedom, enabling the protein to explore its natural conformational ensemble under simulation conditions.
This 1 ns equilibration provides adequate time for structural relaxation while monitoring system stability and thermodynamic properties.
Production¶
The production phase represents the primary data collection period, employing optimized simulation parameters to maximize sampling efficiency while maintaining thermodynamic stability. The protocol employs GPU acceleration for enhanced computational performance during the extended production trajectories.
The GPU-accelerated configuration employs NVIDIA A100 hardware to achieve optimal performance for long-timescale simulations. The single-node configuration minimizes communication overhead while providing substantial computational power for production trajectories.
08_prod_npt¶
The simulation employs the NPT ensemble with a time step of 0.002
ps (2.0
fs) for 50
million MD steps for a total of 100
ns.
Temperature control is achieved through Langevin dynamics with a target temperature of 300
K and a collision frequency of 5.0
ps-1.
Pressure control is implemented using the isotropic position scaling method with a Berendsen barostat, targeting a pressure of 1.01325
bar and a pressure relaxation time of 1.0
ps.
Periodic boundary conditions are used for non-bonded interactions with a 10.0
Å cutoff.
Covalent bonds involving hydrogen are constrained with the SHAKE algorithm.
Coordinates are written every 5000
steps (10
ps) and energy information every 500
steps (1
ps).
YAML
- stage_name: 08_prod_npt
topo_path: "$SLURM_SUBMIT_DIR/inputs/mol.prmtop"
prev_restart_path: "$SLURM_SUBMIT_DIR/../../04-relax/{{ RUN_NAME }}/outputs/07_relax_npt.rst"
ref_coord_path: "$SLURM_SUBMIT_DIR/inputs/mol.inpcrd"
input_kwargs:
irest: 1
ntx: 5
ig: -1
dt: 0.002
nstlim: 50000000
nscm: 500
ntr: 0
ntb: 2
ntf: 2
ntc: 2
cut: 10.0
ntt: 3
temp0: 300.0
gamma_ln: 5.0
ntp: 1
barostat: 2
pres0: 1.01325
mcbarint: 100
comp: 44.6
taup: 1.0
ntxo: 2
ntwr: 5000
ntpr: 500
ntwx: 5000
ioutfm: 1
iwrap: 1
Multiple Replica Strategy¶
The simulation protocol employs three independent replicas to enhance statistical sampling and enable assessment of simulation convergence. Each replica begins from the same minimized and equilibrated structure but employs different random number seeds, ensuring trajectory divergence and independent sampling of conformational space.
The replica-based approach provides several critical advantages for investigating Cu(I) binding mechanisms. First, multiple independent trajectories enable assessment of simulation convergence and identification of reproducible structural features. Second, the enhanced sampling improves statistical precision for calculated ensemble properties such as binding site geometries and hydrogen bonding patterns. Third, the independent trajectories can reveal rare conformational events or binding site fluctuations that might be missed in single replica simulations.