Purpose

This page documents the standard protocol for preparing and running molecular dynamics (MD) simulations of RNA k-turn systems (and other similar small biomolecular solutes) using AMOEBA in Tinker9. It takes a starting PDB structure through system preparation, equilibration, and production MD. Umbrella sampling / pulling is not covered here; see the umbrella sampling / REUS pages for that.

The protocol is implemented by ~/scripts/base_masterscript.sh (copied into each project as masterscript.sh). This page describes what that script does, the steps before it runs, and the reasoning behind each choice.

Overview

Obtain a clean starting PDB
Standardize the PDB via a pdbxyz → xyzpdb round-trip
Solvate + neutralize + build a key file with Tinker-GUI/cli.py
Energy minimization (minimize9)
Heating: 5 temperatures with strong position restraints on the solute (NVT)
Restraint phase-down: gradually relax position restraints at 298 K (NVT)
Production: unrestrained NPT MD

Force field: AMOEBA (amoebabio18.prm). Engine: Tinker9 (dynamic9, minimize9).

1. PDB preparation

Source PDB

Start from a crystal / cryo-EM structure in PDB format. Before running anything:

Decide the protonation states of ionizable residues (for RNA/DNA the standard assumption is deprotonated phosphates; for proteins use PROPKA or similar).
Remove or keep heteroatoms intentionally (crystallographic waters, buffer ions, ligands). pdbxyz will silently drop anything it doesn't recognize.
Decide what to do with missing atoms/residues. If the structure is incomplete, fix it externally first (e.g. with PDBFixer) before the steps below.

pdbxyz → xyzpdb round-trip

The goal of this step is to produce a PDB in the exact format that xyzpdb will emit from then on, so that any future PDB written from an MD trajectory (e.g. for visualization, diffing, or restart) is line-for-line comparable to the starting reference.

pdbxyz  input.pdb  -k tinker.key
xyzpdb  input.xyz  -k tinker.key

Where tinker.key points at the AMOEBA parameter file:

parameters /home/yw24267/programs/tinker/Tinker8/2501/params/amoebabio18.prm

The resulting input.pdb (overwriting the input, so back up the original) is the canonical PDB for this system. Use this as your reference structure from this point on — not the crystal PDB — when comparing to simulation frames.

Why do the round-trip:

The crystal PDB contains headers (HEADER, TITLE, SOURCE, COMPND, REMARK, MASTER, CONECT…) and metadata that xyzpdb does not emit.
Hydrogens added by pdbxyz get standardized names/ordering.
Atom numbering, chain ID handling, and record termination match the Tinker convention after the round-trip.

Without this step, a diff between the input PDB and later xyzpdb output is dominated by format noise rather than real coordinate differences.

2. Solvation and neutralization (Tinker-GUI)

System building is handled by cli.py in Tinker-GUI. The CLI takes a YAML config and produces a minimized-ready .xyz / .key pair plus intermediates in temp/.

config.yaml

Example for an RNA k-turn system with neutralizing Mg²⁺/Na⁺/Cl⁻, no added bulk salt, 12 Å buffer:

tinker_path: /home/yw24267/programs/tinker/Tinker8/2501/bin
amoeba_prm:  /home/yw24267/programs/tinker/Tinker8/2501/params/amoebabio18.prm
output_prefix: minimize
solutes:
  nucleic_acid:
    - myrna.pdb             # use the round-tripped PDB from step 1
solvent:
  name: water
ions:
  neutralizers:
    - Mg2+
    - Na+
    - Cl-
  salts:
    names:
      - Mg2+
      - Na+
      - Cl-
    concentrations:
      - 0.0
      - 0.0
      - 0.0
box:
  type: rectangular          # cubic or rectangular
  buffer: 12.0               # Å between solute and box wall
  # size: [60, 60, 60]       # if set, overrides buffer

Key fields:

output_prefix: minimize — the downstream masterscript expects minimize_final.xyz / minimize_final.key, so keep this as minimize.
neutralizers — the CLI adds whichever of these are needed to zero out the net charge. For RNA (negative backbone) you'll typically need Mg²⁺ and/or Na⁺.
salts.concentrations — set to 0.0 for pure counter-ion neutralization. Use non-zero values (mol/L) to add bulk salt on top of neutralization.
buffer: 12.0 — distance from solute to the nearest box face. 12–15 Å is standard; 15 Å is safer for larger conformational changes.

Running the CLI

python /home/bae969/programs/Tinker-GUI/cli.py -c config.yaml

Outputs:

minimize_final.xyz — solvated, neutralized starting structure
minimize_final.key — matching key file with box dimensions and parameter path
temp/ — intermediate files (solute-only xyz, water box, per-ion-addition snapshots, logs). Do not delete temp/minimize.xyz — the masterscript reads it to determine how many solute atoms to restrain.

Sanity check

Before proceeding, verify the system is sensible:

analyze  minimize_final.xyz  -k minimize_final.key  em

Net charge should be 0 (within rounding). If not, the neutralizer list in config.yaml needs adjustment.
Total dipole should be finite (no NaNs).

analyze  minimize_final.xyz  -k minimize_final.key  ep

Every atom should have multipole and vdW parameters. Missing parameters on any atom means an unrecognized residue / atom name that pdbxyz couldn't map.

3. Minimization

Performed by minimize9 with a gradient tolerance of 0.5 kcal/mol/Å using a key file that includes a position restraint on the solute:

RESTRAIN-POSITION  -1  N  50.0

where N is the number of solute atoms (read by the masterscript from temp/minimize.xyz). The restraint keeps the solute near its starting coordinates while water and ions relax around it. The force constant 50.0 is in kcal/mol/Å².

The exact command run by the masterscript:

minimize9  minimize_final.xyz  -k minimize.key  0.5

Why a restraint during minimization: without it, small clashes between water/ion placement and the solute can pull the solute geometry off the reference structure before dynamics even starts. The restraint also suppresses any collapse of the experimentally-determined conformation during the initial energy drop.

4. Heating (restrained, NVT)

Five sequential NVT runs at increasing temperatures, all with the same key file (heating.key, FC = 50 kcal/mol/Å² on all solute atoms, polar-eps 0.01):

Step	Temperature (K)	Steps	dt (fs)	Total (ps)
Heating_to_10K	10	200 000	2.0	400
Heating_to_75K	75	200 000	2.0	400
Heating_to_150K	150	200 000	2.0	400
Heating_to_225K	225	200 000	2.0	400
Heating_to_298K	298	200 000	2.0	400

Command used for each:

dynamic9  pre_production.xyz  -k heating.key  200000  2.0  10.0  2  <T>  1.0

Argument order: nsteps dt_fs write_interval_ps ensemble(2=NVT) T_K P_atm.

The output trajectory pre_production.arc accumulates across all five heating steps (Tinker9 appends to an existing .arc). Total heating time: 2 ns across 5 temperature rungs.

Why this schedule:

Starting at 10 K avoids the initial velocity kick that would shake apart a minimized-but-cold system.
Going straight to 298 K routinely produces "induced dipole not converging" errors in AMOEBA because of transient bad contacts. Gradual heating in 5 rungs drains that risk.
The strong restraint (FC = 50) keeps the solute essentially frozen while solvent and ions thermalize.

5. Restraint phase-down (NVT, 298 K)

Four sequential NVT runs at 298 K, each stepping the solute position-restraint force constant down:

Step	FC (kcal/mol/Å²)	Steps	Total (ps)
Restraint_10	10.0	200 000	400
Restraint_1	1.0	200 000	400
Restraint_0.1	0.1	200 000	400
Restraint_None	—	200 000	400

Each step uses a separate key file (restraint_10.key, restraint_1.key, restraint_0.1.key, restraint_none.key) that differs only in the RESTRAIN-POSITION line (absent in the last one).

Command:

dynamic9  pre_production.xyz  -k restraint_<FC>.key  200000  2.0  10.0  2  298  1.0

Why step down rather than release at once: an abrupt release at FC = 50 → 0 can cause a sudden "spring release" where the solute snaps to a lower-energy conformation that is different from the experimental reference. A 50 → 10 → 1 → 0.1 → 0 sequence lets the solute relax gradually while solvent/ion distribution tracks it.

After this stage the system is fully equilibrated at 298 K, unrestrained, at the constant-volume box from solvation. 1.6 ns total across the 4 steps.

6. Production (NPT)

The production run is NPT at 298 K, 1 atm, unrestrained, with a tighter polarization convergence:

polar-eps 0.0001 (vs 0.01 during equilibration)
barostat MonteCarlo, thermostat bussi
integrator respa, 2 fs timestep

Starting frame: the masterscript extracts the final frame of the pre_production.arc trajectory (the last xyz_linecount lines) into production.xyz, then launches production from there:

tail -n "$xyz_linecount" pre_production.arc > production.xyz
dynamic9  production.xyz  -k production.key  20000000  2.0  10.0  4  298  1.0

Argument order (NPT): nsteps dt_fs write_interval_ps ensemble(4=NPT) T_K P_atm.

20 000 000 steps × 2 fs = 40 ns production. Frames are written every 10 ps (4 000 frames total). Adjust nsteps for the system / question — 20–100 ns is typical.

Why NPT for production: the box dimensions from the solvation step correspond to a density that is close to but not exactly 1.0 g/cc at 298 K. NPT lets the box settle to the correct density under the AMOEBA potential. Skipping this and running NVT with the as-built box can leave the system slightly compressed or stretched, which biases structural observables.

Why tighter polar-eps: 0.01 is fine for equilibration (faster, noise is averaged out over the short run). For production the induced-dipole error should be as small as practical because it enters every observable; 1e-4 is the routine choice and 1e-5 is worth it if energetics are the main observable.

7. Running the whole pipeline

Copy ~/scripts/base_masterscript.sh into the project directory as masterscript.sh. It assumes:

minimize_final.xyz and minimize_final.key are present in the current directory (outputs of step 2).
temp/minimize.xyz is present (so N = solute atom count can be read).
Tinker9 binaries are at /home/yw24267/programs/tinker/Tinker9/latest/build (edit TINKER_BIN if yours differ).
The correct GPU is selected via CUDA_VISIBLE_DEVICES.

Typical launch on a GPU node:

ssh nodeXXX
cd /path/to/project
export CUDA_VISIBLE_DEVICES=0
nohup bash masterscript.sh > masterscript.log 2>&1 &

What the script generates on the fly (overwriting each run):

heating.key, minimize.key (FC = 50)
restraint_10.key, restraint_1.key, restraint_0.1.key
restraint_none.key, production.key

All derived from minimize_final.key — it reads the a-axis / b-axis / c-axis / pme-grid lines from there and appends them, then adds or omits the RESTRAIN-POSITION line and sets the appropriate polar-eps.

Summary of wall-clock and file sizes

Stage	MD time	Runtime (rough, single GPU)	Output
Minimization	—	minutes	minimize_final.xyz_2
Heating (5 × 400 ps)	2.0 ns	~few hours	pre_production.arc (appended)
Restraint phase-down	1.6 ns	~few hours	pre_production.arc (appended)
Production (20 M steps)	40 ns	~1–2 days	production.arc

Runtime scales with system size and GPU generation; the numbers above are order-of-magnitude for a ~700-atom solute in a ~60 Å water box on a modern GPU.

Notes and common pitfalls

Induced dipole not converging during heating almost always means the starting structure has a bad contact. Re-check minimization went to a reasonable gradient, and that the solvation step didn't place a water atop a solute atom. You can temporarily add USOLVE-CUTOFF 0.0 to the heating key if convergence is marginal; remove for production.
Net charge non-zero after solvation means the neutralizers list in config.yaml didn't cover the real net charge. Re-run with the right ion set.
RESTRAIN-POSITION force constant is in kcal/mol/Å², not kcal/mol. The masterscript embeds this unit convention; do not change the numbers without understanding the scale.
Do not delete temp/ until after the masterscript has read temp/minimize.xyz. If you've deleted it, re-running Tinker-GUI/cli.py regenerates it.
Restarting a crashed run: Tinker9 resumes from .dyn if present. To restart from scratch, remove both .dyn and .arc (otherwise the new trajectory appends to the old one).
Stopping a running MD cleanly: touch pre_production.end (or production.end). Tinker9 checks for this file each frame and stops gracefully, leaving a valid .dyn to resume from.

Tutorial:kturn tut

Contents

Purpose

Overview

1. PDB preparation

Source PDB

pdbxyz → xyzpdb round-trip

2. Solvation and neutralization (Tinker-GUI)

config.yaml

Running the CLI

Sanity check

3. Minimization

4. Heating (restrained, NVT)

5. Restraint phase-down (NVT, 298 K)

6. Production (NPT)

7. Running the whole pipeline

Summary of wall-clock and file sizes

Notes and common pitfalls

See also

Navigation menu

Tutorial:kturn tut

Purpose

Overview

1. PDB preparation

Source PDB

pdbxyz → xyzpdb round-trip

2. Solvation and neutralization (Tinker-GUI)

config.yaml

Running the CLI

Sanity check

3. Minimization

4. Heating (restrained, NVT)

5. Restraint phase-down (NVT, 298 K)

6. Production (NPT)

7. Running the whole pipeline

Summary of wall-clock and file sizes

Notes and common pitfalls

See also

Navigation menu

Search