Usage#
If you are not yet familiar with GROMACS, please follow some tutorials of the official GROMACS before starting the SAXS/WAXS calculations. You should be familiar with normal MD simulations with GROMACS and learn about the relevant files (top, itp, mdp, ndx, tpr, etc.). You will find several GROMACS tutorials with your favourite search engine, for instance the excellent tutorials by the Justin Lemkuhl lab.
An introduction into explicit-solvent SAXS calculations is available as preprint here. An introduction to the SAXS-driven MD simulations is available here.
Modified and added GROMACS modules#
With SAS (small-angle scattering), we refer in the following to SAXS, WAXS, or SANS curves.
gmx mdrun -rerun: SAS curves are computed from a given trajectory with the mdrun module using the rerun functionality [CH14]. The SAS calculations are activated by mdp options (see below).
gmx mdrun: mdrun is also used for running a SAS-driven simulations, that is, simulations that are coupled to one or multiple (experimental) target curves [CH15, CSS+19]. Such simulations are used to refine a structure against experimental data. The (experimental) target SAS curve(s) are provided with option gmx mdrun -is.
gmx genscatt: The new genscatt module writes itp files that contain definitions for atomic form factors (so-called Cromer-Mann parameters), used to compute X-ray scattering amplitudes. Optionally, genscatt adds neutron scattering lengths definitions to the itp file for computing SANS curves. Check gmx genscatt -h for details. The written itp file, e.g. scatter.itp, must be included into the [ moleculetype ] definition, similar to the posre.itp for defining position restraints.
Envelope around a GB3 domain. Atoms inside the envelope contribute to the SAS predictions.#
gmx genenv: The new genenv modules generates the envelope around the biomolecule. All solvent atoms within the envelope contribute to the SWAXS predictions. To fully account for the hydration layer, the envelope should typically be constructed at a distance of 6-7 Angstrom from the biomolecule. Check gmx genenv -h for details. Take a look at the envelope within PyMol with run envelope.py
gmx mdrun input options#
The following gmx mdrun options are relavant for SAS predictions and the SAS-driven MD simulations:
-ow (waxs). Top name for SAS-related output files. Output will be written to files such as waxs_final.xvg, waxs_contrib.xvg, waxs.log etc.
-sw (water.tpr). Run input file (tpr) of the pure-solvent simulation used to compute the buffer subtraction. The tpr file must be generated with grompp from this modified GROMACS version. To make sure that the water.tpr file contains atomic form factor definitions (Cromer-Mann parameters) of the solvent atoms, the tpr must be created with an mdp file with options waxs-solvent=System, while leaving waxs-solute empty). Make sure that the pure-solvent system is large enough such that the envelope fits into the compact unit cell. In doubt, just make the pure-solvent box a little larger than the solute box.
-fw (water.xtc). XTC trajectory file of pure-solvent simulation. The number of frames used from the XTC is controlled with mdp option waxs-nfrsolvent. To ensure that the buffer matches between solute and pure-solvent systems, the solvent should have been simulated with the same force field and mdp options as the solute.
-is (target.xvg). The experimental target curve for SAS-driven MD simulations. The file target.xvg must have 3 columns (q, Iexp, \({\sigma}\)). If you don’t have errors \({\sigma}\), you can simulate realistic errors following Sedlak et al. [SBL17] and add them to the target.xvg. For a quick test, you could also use an error of 1% of Iexp.
gmx mdrun output#
gmx mdrun writes the following SWAXS-related output:
waxs_final.xvg. The most important output: the calculated SWAXS or SANS curve Icalc(q). If you requested multiple SAXS and SANS curves, these files are numbered as waxs_final_0.xvg, waxs_final_1.xvg etc.
prot+solvlayer_0.pdb. PDB structure of the solute and all solvent atoms inside the envelope, which contribute to the SAS curve, taken from the first frame. See the image on the right. Take a look in PyMol together with the envelope file envelope.py. More such frames are only written if the verbosity level is increased with export GMX_WAXS_VERBOSE=2 or higher.
excludedvolume_0.pdb. PDB structure of the excluded solvent of the first frame, as defined by the envelope. Again, more frames are written with a larger verbosity level.
waxs_pot.xvg. SAXS-derived potential energy vs. simulation timee, shown for each q-point. This shows which
q-range mainly drives a conformational transition in a SAXS-driven simulation.
waxs.log. Log file, most relevant are the statistics (including the radius of gyration) at the end.
waxs_contrib.xvg. The mathematical contributions to \({I_\mathrm{calc}(q)}\), see Ref. [CH14] for definitions. These are not relevant for most users.
waxs_spectra.xvg. \({I_\mathrm{calc}(q)}\) versus simulation time, one curve every waxs-nstlog. This shows how the SAS curve converges with increasing number of simulation frames.
waxs_stddevs.xvg. Statistical and systematic errors (\({\sigma_\mathrm{exp}}\), \({\sigma_\mathrm{calc}}\), \({\sigma_\mathrm{buf}}\)).
waxs_averageA.xvg, waxs_averageB.xvg, waxs_averageD.xvg.
Average scattering amplitude of the Protein with solvation layer \({\langle A_i(\mathbf{q})\rangle}\),
the excluded solvent \({\langle B_i(\mathbf{q})\rangle}\),
and of the buffer-subtracted scattering amplitude \({\langle D(\mathbf{q})\rangle}\).
See Ref. [CH14] for a detailed explanation. These can, for instance, be used to compute anisotropic scattering patters. Not relevant for normal users.
waxs_density.dat, waxs_density.dat.cube. Electron density within the envelope, written only with environment variable GMX_WAXS_GRID_DENSITY=1, see below. Cube format can be visualized by, e.g., Pymol.
mdp options#
The following mdp options control the SAS predictions and the SAS-driven MD simulations.
define. Keep define = -DSCATTER in the mdp file for all SAS calculations since the scattering definitions in the topology and force field files are defined in blocks #ifdef SCATTER ... #endif. If you want to use your topology with official GROMACS, simply remove -DSCATTER from your mdp file.
scatt-coupl (xray/neutron). For computing a SWAXS curve, or for coupling to a SWAXS curve use xray. For computing SANS curves, or for coupling to a SANS curve, use neutron. Multiple options are supported, e.g. xray neutron neutron for computing the SAXS curve and two SANS curves at two different D2O concentrations. For each neutron keyword, one value at waxs-deuter-conc specifies the D2O concentration.
waxs-solute. Group of solute’s physical atoms used for SWAXS/SANS calculations, e.g. Protein. If your solute contains virtual sites, use Prot-Masses. If your solute is not a default index group generated by GROMACS, (such as Protein with a heme group, Protein/DNA complex, Protein with a bound ion, membrane protein in detergent belt), prepare an appropriate index group with the GROMACS module gmx make_ndx or gmx select.
waxs-solvent. Group of solvent atoms, e.g. Water or Water_and_ions. If you have salt in your solvent, make sure to include it in waxs-solvent
waxs-rotfit. Group used for rotational fit of the solute into the shaped envelope. This group must correspond to the rotational fit group you specified with gmx genenv. A natural choice would be C-alpha.
waxs-pbcatom (-1). The atom ID used to make the solute whole, similar to the PBC atom used for pull groups. Special allowed values are:
-1: use number-wise center-2: The atomic distances between numberwise atomic neighbors are used to make the solute whole. Hence,-2typically works fine if your solute contains only a single chain, but-2does not work with a multi-domain protein.
The PBC atom should be reasonably near the geometric center of the solute. gmx genenv suggests a good PBC atom. The atom ID is the global atom ID, not the solute-internal ID. If you use an unsuitable PBC atom, you will get an error message, but it will not influence the result of the calculation.
waxs-tau (ps).
Memory time \({\tau}\) for SWAXS calculations in picoseconds.
-1: uniformly-weighted average. Use this when computing a SWAXS curve from a given simulations withgmx mdrun -rerun.For a SAS-driven simulation, use a value \({\tau}\)>0 - this way an on-the-fly running average is computed with weights (or a memory kernel) that decay exponentially into the past. This averageing on-the-fly ensures that the SAS curve is reasonably converged before it is compared with the target curve provide with
gmx mdrun -is. A reasonable value for SAXS-driven simulations is between 50ps and 500ps.
waxs-t-target (ps).
During SAXS-driven simulations, turn on the SAXS-derived restraints
gradually within this time (e.g. 10000 ps). This avoids strong SAXS-derived forces at the beginnig of the simulation when the solute is in a conformation that is not compatible with the given target SAS curve. With waxs-t-target, you allow the
solute carry out the transition more slowly.
waxs-potential (log/linear).
linear: Recommended! Apply a harmonic restraint on the intensity I(q). Hence, the SAS-derived potential follows \({V(\mathbf{r}) \propto \sum_i (I_\mathrm{calc}(q_i) - \tilde{I}_\mathrm{exp}(q_i))^2/\sigma_i^2}\). With this coupling scheme, it you should usewaxs-weights=exp+calc+solvdens, for the following reason: The optionsolvdensindicates that our error model accounts for systematic errors due to an uncertainty from the buffer subtraction, thereby reducing the weights of the SAS-derived forces at very small q, where the true uncertainty is dominated by systematic errors and not by statistical errors. Otherwise, you may get spruriously high forces from the very small q-region and your simulation may crash.log: Apply a harmonic restraint on the logarithm of the intensity I(q). The SAS-derived potential follows \({V(\mathbf{r}) \propto \sum_i (\ln I_\mathrm{calc}(q_i) - \ln \tilde{I}_\mathrm{exp}(q_i))^2/(\sigma(\ln I_\mathrm{calc}(q_i))^2}\). Here, the uncertainty of the logarithm of the intensity is \({\sigma(\ln I_\mathrm{calc}(q_i)) = \sigma_i/I_\mathrm{calc}(q_i)}\). With this coupling scheme, typically usewaxs-weights=uniform. The restraint on the logarithm of I(q) leads to roughly uniform weights along the q region, and a reasonable balance between contributions from smaller and wider q regions.
waxs-weights (uniform, exp, exp+calc, exp+solvdens, or exp+calc+solvdens).
Specify how the weights of q-points durig SAXS-driven MD are computed.
uniform: equal weights on all q-points, no uncertainties \({\sigma_i}\) are applied. Only useful withwaxs-potential=logexp+calc+solvdens: Taking uncertainty from the experimental curve \({\sigma_\mathrm{exp}(q_i)}\), from the caluculated curve \({\sigma_\mathrm{c}(q_i)}\), and systematic error from the buffer subtraction \({\sigma_\mathrm{buf}(q_i)}\) into account. The total uncertainty is computed via \({\sigma^2(q_i) = \sigma_\mathrm{exp}^2(q_i) + \sigma_\mathrm{c}^2(q_i) +\sigma_\mathrm{buf}^2(q_i) }\). Seewaxs-solvdens-uncertfor controlling \({\sigma_\mathrm{buf}(q_i)}\). This is recommended withwaxs-potential=linear.exp,calc,exp+solvdens: Use only certain sources of errors. These options are not recommended.
waxs-fc.
Force constant for SAXS-driven MD. With waxs-potential=linear, a good value for waxs-fc is around 1. With waxs-potential=log, a good value for waxs-fc is around 300. To find a good value for
waxs-fc, we recommend to carry out multiple SAXS-driven MD imulations with different
force constants, and to observe the behaviour of the biomolecule. Is the transition
plausible? Is the final SAS curve in agreement with the target curve? If waxs-fc too small, no transition may occur and the computed SAS curve does not match the target curve. If waxs-fc is too large, your protein might start to unfold or the simulation may become unstable.
Bayesian interpretation: With waxs-fc=1 and waxs-potential=linear, the ensemlbe generated by your SAS-driven MD simulation can be interpreted as the posterior distribution of a Bayesian inference problem, see [SH17]. Hence, the width of your ensemle provides the uncertainty of the refined structure in the light of the SAS data and the force field. In the context of NMR refinement, this method was referred to as Inferential Structure Determination (ISD) [RHN05].
Hint: In case you need a larger waxs-fc to drive the conformational transition, you may switch to waxs-fc=1 in a follow-up simulation, thereby allowing the (elegant) Bayesian interpretation.
waxs-nstcalc (steps).
Frequency of updating the SAXS curve (in number of MD steps). For computing a
SAXS curve from a given simulations with gmx mdrun -rerun,
simply use waxs-nstcalc=1. In a SWAXS-driven simulation,
choose waxs-nstcalc such that the SAXS curve is computed
from statistically independent solvent configurations, for instance every 0.5ps. With integration time steps of
2 fs or 4 fs, this would be given with waxs-nstcalc=250 or waxs-nstcalc=125, respectively.
waxs-nfrsolvent.
Number of pure-solvent frames used to compute the scattering of the
excluded solvent, denoted B(q) in [CH14]. At the beginning of a simulation or a rerun,
the scattering from one solvent frame B(q) is computed per scattering from each solute frame (denoted A(q)). After waxs-nfrsolvent SAXS steps, B(q) is not updated any more. For many applications, waxs-nfrsolvent=1000 is save choice.
For system with a small electron density contrast between solute and solvent (such as intrinsically disordered proteins),
you may need to use larger values to ensure that the SAXS curve is converged at small q values. For large systems with a
large contrast a much smaller values may suffice. In doubt, compute the SAXS from indpendent solute and pure-solvent frames multiple times
and compare the results.
waxs-nstlog (steps):
Frequency of output to waxs_spectra.xvg and to waxs_stddevs.xvg.
waxs-nq
Number of q-points in calculated SAS curve \({I_\mathrm{calc}(q_i)}\).
SAS calculation with
gmx mdrun -rerunfrom a given trajectory :waxs-nq=101is reasonable.SAXS-driven MD simulations: Since (theoretical) SAS curves are smooth, there is no point using many more q-points compared to the number of Shannon-Nyqvist channels \({N_{sn}}\). Hence, we recommend computing \({N_{sn} = (q_\mathrm{max} - q_\mathrm{min})\cdot D_\mathrm{sol}/\pi}\), where \({q_\mathrm{max}}\) and \({q_\mathrm{min}}\) denote the maximum and minimum q-values (see
waxs-startqandwaxs-endq), and \({D_\mathrm{sol}}\) is the maximum diameter of the solute. Then use a somewhat larger values such as \({2\times N_{sn}}\) (e.g.waxs-nq=20if Nsn=10).
Important: Raw exeprimental \({I_{raw}(q)}\) curves are often heavily oversampled and exhibit large statisical errors. In other words, there are many experiental data points per Shannon-Nyqvist channel. When doing SAXS-driven MD simulations by coupling to \({\approx 2\times N_{sn}}\) q-points, it is critical to first smooth the experimental curve \({I_{raw}(q)}\) and to estimate the uncertainty of the \({I_{exp}(q)}\). We provide a shell script smooth-saxs-curve.sh that
uses the ATSAS modules DATGNOM and DATREGRID for smoothing and re-binning the SAS curve. The target curve target.xvg provided with gmx mdrun -is must contain the smoothed and re-binned experimental curve.
waxs-startq:
starting momentum transfer q (inverse nanometer nm-1, not inverse Angstrom)
waxs-endq: largest q (nm-1)
waxs-nsphere:
Number of q-vectors used per absolute value of q for computing
the orientational average (spherical quadrature), denoted \({J}\) in [CH14].
A appropriate waxs-nsphere value strongly depends on the solute diameter \({D_\mathrm{sol}}\)
and the maximum \({q_\mathrm{max}}\); we recommend using at least \({J > 0.2\cdot(q_\mathrm{max}D_\mathrm{sol})^2}\), where
\({q_\mathrm{max}}\) is the maximum momentum transfer and \({D_\mathrm{sol}}\) the diameter of the solute.
Example: For a solute with diameter 5 nm and a SAS curve up to q = 10nm-1, use at least waxs-nsphere=500. In doubt, try a few different values.
waxs-Iexp-fit (no/scale-and-offset/scale):
During SAS-driven MD, fit the experimental SAS curve \({I_\mathrm{exp}(q)}\) against the calculated curve
\({I_\mathrm{calc}(q)}\) before computing the SAS-derived potential and forces. This is done to remove
the overall scale, which is essentially arbitrary, and, optionally, a constant offset to absorb
some uncertainty due to the buffer subtraction.
scale-and-offset: Fit experimental curve as \({I_\mathrm{fit}(q) = fI_\mathrm{exp}(q)} + c\). This is recommended since uncertainties in the buffer subtraction are quite common.scale: Fit experimental curve as \({I_\mathrm{fit}(q) = fI_\mathrm{exp}(q)}\)no: Do not fit the experimental curve.
waxs-solvdens:
Experimental electron density of the solvent (e/nm3).
Default: 334. The value is used for the solvent density correction. Such correction is needed because many
popular water model do not reproduce the experiental electron denisty. For simulations at different temperatures
or for simulations containing heavy ions, waxs-solvdens must be adapted to match the experimental density.
waxs-solvdens-uncert:
Relative uncertainty of the experimental buffer density,
used to estimate a systematic error \({\sigma_\mathrm{buf}}\) due to the buffer subtraction. A reasonable
value is between 0.001 (0.1%) and 0.005 (0.5%). Only applied if waxs-weights contains the option solvdens.
Larger values of waxs-solvdens-uncert reduce the weights at very small q-region. Hence, if your experimental
target curve is compromised at small q due to systematic errors from aggregation or particle repulsion, use a slightly
larger value such as 0.01. In doubt, have a look at waxs_pot.xvg to see how larger waxs-solvdens-uncert lead to
smaller SAXS-derived
potential energies at small q.
waxs-solvdens-uncert-bayesian (no/yes):
Sample the systematic error of the solvent density
in a Bayesian manner, simultaneously with the refined structure (see Shevchuk and Hub, PLoS Comput Biol, 2017).
For the uncertainty of the
solvent density, a Gaussian prior is assumed with a width defined by waxs-solvdens-uncert. The sampled
uncertainties are written into the file waxs_solvDensUncert.xvg. A histogram over the values in waxs_solvDensUncert.xvg
will give you the posterior distribution of the systematic error. The larger the systematic errors in the experimental curve at
low q, the more this posterior will shift to larger values.
waxs-correct-buffer (no/yes).
Specifies the convention of buffer subtraction.
no: \({I(q) = I_\mathrm{sample}(q) - I_\mathrm{buffer}(q)}\) (default, most common buffer subtraction scheme in the SAS community)yes: \({I(q) = I_\mathrm{sample}(q) - (1-v)I_\mathrm{buffer}(q)}\) where \({v}\) denotes the volume fraction of the solute in the solution. The correction to the SAXS curve withwaxs-correct-buffer=yesis computed following [KofingerH13].
waxs-deuter-conc:
Specify the deuterium concentration of the buffer, one value for each keyword neutron
in the mdp option scatt-coupl.
Example: with scatt-coupl = xray neutron neutron neutron you may use waxs-deuter-conc = 0 0.42 1 for computing the SANS curves at D2O concentrations of 0%, 42% and 100%.
waxs-scale-i0 (no/yes):
During SAS-driven MD, correct the overall electron density contrast such that \({I_\mathrm{calc}(q=0)}\) matches a given value.
The given value is either taken from the target curve target.xvg, or from the value of the environment variable GMX_WAXS_I0 (if defined, example: export GMX_WAXS_I0=3.56e6). Here, the target curve or the value of GMX_WAXS_I0 must be given in units \({e^2}\) as used by GROMACS-SWAXS.
Matching \({I_\mathrm{calc}(q=0)}\) is achieved by adding a small uniform electron density to the excluded solvent. This feature is useful
if the forward scattering converges slowly, as common for intrinsically disordered proteins (IDPs) that exhibit only a small contrast relative to the buffer [HH19]. For IDPs, \({I_\mathrm{calc}(q=0)}\) does not converge within the memory time waxs-tau and, hence, could not be compared with the target curve in a meaningful way. With, waxs-scale-i0=yes, in contrast, the convergence of \({I_\mathrm{calc}(q=0)}\) is greatly accellerated, allowing for SAXS-driven MD simulations [HH19].
waxs-warnlay (nm).
mdrun writes a warning if the solute is getting
closer than the distance waxs-warnlay to the surface of the envelope (deprecated).
waxs-anisotropic: deprecated, removed in GROMACS-SWAXS 2021.4-0.2
waxs-energy: deprecated, removed in GROMACS-SWAXS 2021.4-0.3
waxs-ensemble-type, waxs-ensemble-nstates, waxs-ensemble-init-w, waxs-ensemble-fc : Options for ensemble simulations. Documentation coming soon.
Environment variables#
The behavior of gmx mdrun can be modfied by the following of environment variables:
Mandatory environment variables for gmx mdrun#
GMX_ENVELOPE_FILE Path to the envelope file witten by gmx genenv (such as envelope.dat)
GMX_WAXS_FIT_REFFILE Path to the envelope reference coordinate file written by gmx genenv (such as envelope-ref.gro), used to fit the solute into the envelope
Other useful environment variables#
Boolean variables are turned on by setting them to 1. Example: export GMX_WAXS_GRID_DENSITY=1
GMX_WAXS_BEGIN, GMX_WAXS_END (ps). Start and end time in picoseconds for computing SWAXS/SANS curves from given trajectory with gmx mdrun -rerun traj_comp.xtc, same as options -b and -e in many GROMACS tools.
GMX_WAXS_GRID_DENSITY. If set to 1, write average electron denisty withing the envelope to waxs_density.dat and in Cube format to waxs_density.dat.cube. The Cube format can, for instance, be visualized by Pymol.
GMX_WAXS_GRID_DENSITY_MODE (0,1,2). Control, whether density of solute, solvent, or both are written to the grid density file.
0: density of solute and solvent (default)1: density of solute only2: density of solvent only
Example: To write the electon density of the solvent to a Cube file, set: export GMX_WAXS_GRID_DENSITY=1; export GMX_WAXS_GRID_DENSITY_MODE=2
WAXS_ENVELOPE_GRID_DENSITY_SPACING. Spacing of grid density. Default: 1 Angström.
More specialized environment variables#
Boolean variables are turned on by setting them to 1, true, or yes.
GMX_WAXS_PRINTFORCES. Write SAXS-derived forces into a trr file waxs_forces.trr.
GMX_WAXS_VERBOSE. Set verbosity level (0, 1, 2, or 3)
GMX_WAXS_SOLUTE_IS_WHOLE. Do not try to make the solute whole using the waxs-pbc atom. Assume that the solute is already whole in the trajectory provided with mdrun -rerun.
GMX_WAXS_HAVE_FITTED_TRAJ. Do not fit the solute into the envelope using reference file (envelope-ref.gro), instead assume that the solute is already fitting into the envelope. Together with GMX_WAXS_SOLUTE_IS_WHOLE, this can be useful if it is difficult to make the molecule whole and to fit it. For example, if your solute is a virus capsid (sphere with a large hole), there is no good choice for a waxs-pbc atom.
GMX_WAXS_PRINT_NINDEP. Write the estimated number of independent q-vectors per absolute value of q to waxs_nindep.xvg. These are used to compute the statistical errors of Icalc(q) when doing the spherical quadratures (spherical average) over q-vectors.
References#
Po-chia Chen and Jochen S. Hub. Validating solution ensembles from molecular dynamics simulation by wide-angle X-ray scattering data. Biophys. J., 107:435–447, 2014.
Po-chia Chen and Jochen S. Hub. Interpretation of solution X-ray scattering by explicit-solvent molecular dynamics. Biophys. J., 108:2573–2584, 2015.
Po-chia Chen, Roman Shevchuk, Felix M. Strnad, Charlotte Lorenz, Lukas Karge, Ralph Gilles, Andreas M. Stadler, Janosch Hennig, and Jochen S. Hub. Combined Small-Angle X-ray and Neutron Scattering Restraints in Molecular Dynamics Simulations. Journal of Chemical Theory and Computation, 15(8):4687–4698, August 2019. doi:10.1021/acs.jctc.9b00292.
Markus R. Hermann and Jochen S. Hub. SAXS-Restrained Ensemble Simulations of Intrinsically Disordered Proteins with Commitment to the Principle of Maximum Entropy. Journal of Chemical Theory and Computation, 15(9):5103–5115, September 2019. doi:10.1021/acs.jctc.9b00338.
Jürgen Köfinger and Gerhard Hummer. Atomic-resolution structural information from scattering experiments on macromolecules in solution. Phys. Rev. E., 87:052712, 2013.
Wolfgang Rieping, Michael Habeck, and Michael Nilges. Inferential structure determination. Science, 309(5732):303–306, 2005.
Steffen M. Sedlak, Linda K. Bruetzel, and Jan Lipfert. Quantitative evaluation of statistical errors in small-angle X-ray scattering measurements. J Appl Crystallogr, 50(2):621–630, April 2017. doi:10.1107/S1600576717003077.
Roman Shevchuk and Jochen S. Hub. Bayesian refinement of protein structures and ensembles against SAXS data using molecular dynamics. PLoS Comput. Biol., 13:e1005800, 2017.