<li><A HREF="#general"><b>General remarks</b></A>
<p> </p>
<li><A HREF="#pp"><b>preprocessing</b></A> (include, define)
-<li><A HREF="#run"><b>run control</b></A> (integrator, tinit, dt, nsteps, init_step, comm_mode, nstcomm, comm_grps)
-<li><A HREF="#ld"><b>langevin dynamics</b></A> (bd_fric, ld_seed)
+<li><A HREF="#run"><b>run control</b></A> (integrator, tinit, dt, nsteps, init-step, comm-mode, nstcomm, comm-grps)
+<li><A HREF="#ld"><b>langevin dynamics</b></A> (bd-fric, ld-seed)
<li><A HREF="#em"><b>energy minimization</b></A> (emtol, emstep, nstcgsteep)
<li><a HREF="#xmdrun"><b>shell molecular dynamics</b></a>(emtol,niter,fcstep)
<li><a HREF="#tpi"><b>test particle insertion</b></a>(rtpi)
-<li><A HREF="#out"><b>output control</b></A> (nstxout, nstvout, nstfout, nstlog, nstcalcenergy, nstenergy, nstxtcout, xtc_precision, xtc_grps, energygrps)
-<li><A HREF="#nl"><b>neighbor searching</b></A> (nstlist, ns_type, pbc, periodic_molecules, rlist, rlistlong)
-<li><A HREF="#el"><b>electrostatics</b></A> (coulombtype, rcoulomb_switch, rcoulomb, epsilon_r, epsilon_rf)
-<li><A HREF="#vdw"><b>VdW</b></A> (vdwtype, rvdw_switch, rvdw, DispCorr)
-<li><A HREF="#table"><b>tables</b></A> (table-extension, energygrp_table)
-<li><A HREF="#ewald"><b>Ewald</b></A> (fourierspacing, fourier_nx, fourier_ny, fourier_nz, pme_order, ewald_rtol, ewald_geometry, epsilon_surface, optimize_fft)
-<li><A HREF="#tc"><b>Temperature coupling</b></A> (tcoupl, nsttcouple, tc_grps, tau_t, ref_t)
+<li><A HREF="#out"><b>output control</b></A> (nstxout, nstvout, nstfout, nstlog, nstcalcenergy, nstenergy, nstxtcout, xtc-precision, xtc-grps, energygrps)
+<li><A HREF="#nl"><b>neighbor searching</b></A> (nstlist, ns-type, pbc, periodic-molecules, rlist, rlistlong)
+<li><A HREF="#el"><b>electrostatics</b></A> (coulombtype, rcoulomb-switch, rcoulomb, epsilon-r, epsilon-rf)
+<li><A HREF="#vdw"><b>VdW</b></A> (vdwtype, rvdw-switch, rvdw, DispCorr)
+<li><A HREF="#table"><b>tables</b></A> (table-extension, energygrp-table)
+<li><A HREF="#ewald"><b>Ewald</b></A> (fourierspacing, fourier-nx, fourier-ny, fourier-nz, pme-order, ewald-rtol, ewald-geometry, epsilon-surface, optimize-fft)
+<li><A HREF="#tc"><b>Temperature coupling</b></A> (tcoupl, nsttcouple, tc-grps, tau-t, ref-t)
<li><A HREF="#pc"><b>Pressure coupling</b></A> (pcoupl, pcoupltype,
- nstpcouple, tau_p, compressibility, ref_p, refcoord_scaling)
-<li><A HREF="#sa"><b>simulated annealing</b></A> (annealing, annealing_npoints, annealing_time, annealing_temp)
-<li><A HREF="#vel"><b>velocity generation</b></A> (gen_vel, gen_temp, gen_seed)
-<li><A HREF="#bond"><b>bonds</b></A> (constraints, constraint_algorithm, continuation, shake_tol, lincs_order, lincs_iter, lincs_warnangle, morse)
-<li><A HREF="#egexcl"><b>Energy group exclusions</b></A> (energygrp_excl)
-<li><A HREF="#walls"><b>Walls</b></A> (nwall, wall_type, wall_r_linpot, wall_atomtype,
-wall_density, wall_ewald_zfac)
+ nstpcouple, tau-p, compressibility, ref-p, refcoord-scaling)
+<li><A HREF="#sa"><b>simulated annealing</b></A> (annealing, annealing-npoints, annealing-time, annealing-temp)
+<li><A HREF="#vel"><b>velocity generation</b></A> (gen-vel, gen-temp, gen-seed)
+<li><A HREF="#bond"><b>bonds</b></A> (constraints, constraint-algorithm, continuation, shake-tol, lincs-order, lincs-iter, lincs-warnangle, morse)
+<li><A HREF="#egexcl"><b>Energy group exclusions</b></A> (energygrp-excl)
+<li><A HREF="#walls"><b>Walls</b></A> (nwall, wall-type, wall-r-linpot, wall-atomtype,
+wall-density, wall-ewald-zfac)
<li><A HREF="#pull"><b>COM pulling</b></A> (pull, ...)
-<li><A HREF="#nmr"><b>NMR refinement</b></A> (disre, disre_weighting, disre_mixed, disre_fc, disre_tau, nstdisreout, orire, orire_fc, orire_tau, orire_fitgrp, nstorireout)
-<li><A HREF="#free"><b>Free Energy calculations</b></A> (free_energy, init_lambda, delta_lambda, sc_alpha, sc_power, sc_sigma, couple-moltype, couple-lambda0, couple-lambda1, couple-intramol)
-<li><A HREF="#neq"><b>Non-equilibrium MD</b></A> (acc_grps, accelerate, freezegrps, freezedim, cos_acceleration, deform)
-<li><A HREF="#ef"><b>Electric fields</b></A> (E_x, E_xt, E_y, E_yt, E_z, E_zt )
+<li><A HREF="#nmr"><b>NMR refinement</b></A> (disre, disre-weighting, disre-mixed, disre-fc, disre-tau, nstdisreout, orire, orire-fc, orire-tau, orire-fitgrp, nstorireout)
+<li><A HREF="#free"><b>Free energy calculations</b></A> (free-energy, nstfep, nstdgdl, dhdl-print-energy, init-lambda, delta-lambda, fep-lambdas, coul-lambdas, vdw-lambdas, bonded-lambdas, restraint-lambdas, mass-lambdas, sc-alpha, sc-coul, sc-power, sc-r-power, sc-sigma, couple-moltype, couple-lambda0, couple-lambda1, couple-intramol)
+<li><A HREF="#expanded"><b>Expanded ensemble simulation</b></A> (lmc-stats, lmc-mc-move, lmc-seed, lmc-gibbsdelta, mc-temperature, nst-transition-matrix,init-lambda-weights,initial-wl-delta,wl-scale,wl-ratio,symmetrized-transition-matrix,lmc-forced-nstart,weight-c-range,mininum-var-min,lmc-weights-equil,weight-equil-wl-delta,weight-equil-number-all-lambda,weight-equil-number-steps,weight-equil-number-samples,weight-equil-count-ratio,simulated-tempering,simulated-tempering-scaling,sim-temp-low,sim-temp-high)
+<li><A HREF="#neq"><b>Non-equilibrium MD</b></A> (acc-grps, accelerate, freezegrps, freezedim, cos-acceleration, deform)
+<li><A HREF="#ef"><b>Electric fields</b></A> (E-x, E-xt, E-y, E-yt, E-z, E-zt )
<li><A HREF="#qmmm"><b>Mixed quantum/classical dynamics</b></A> (QMMM, QMMM-grps, QMMMscheme, QMmethod, QMbasis, QMcharge, Qmmult, CASorbitals, CASelectrons, SH)
-<li><A HREF="#gbsa"><b>Implicit solvent</b></A> (implicit_solvent, gb_algorithm, nstgbradii, rgbradii, gb_epsilon_solvent, gb_saltconc, gb_obc_alpha, gb_obc_beta, gb_obc_gamma, gb_dielectric_offset, sa_algorithm, sa_surface_tension)
-<li><A HREF="#user"><b>User defined thingies</b></A> (user1_grps, user2_grps, userint1, userint2, userint3, userint4, userreal1, userreal2, userreal3, userreal4)
+<li><A HREF="#gbsa"><b>Implicit solvent</b></A> (implicit-solvent, gb-algorithm, nstgbradii, rgbradii, gb-epsilon-solvent, gb-saltconc, gb-obc-alpha, gb-obc-beta, gb-obc-gamma, gb-dielectric-offset, sa-algorithm, sa-surface-tension)
+<li><A HREF="#adress"><b>AdResS settings</b></A> (adress, adress_type, adress_const_wf, adress_ex_width, adress_hy_width, adress_ex_forcecap, adress_interface_correction, adress_site, adress_reference_coords, adress_tf_grp_names, adress_cg_grp_names)
+<li><A HREF="#user"><b>User defined thingies</b></A> (user1-grps, user2-grps, userint1, userint2, userint3, userint4, userreal1, userreal2, userreal3, userreal4)
<li><A HREF="#idx"><b>Index</b></A>
</ul>
</P>
Depending on the computational cost of the force calculation,
this can take a significant part of the simulation time.
The temperature for one or more groups of atoms
-(<b><A HREF="#tc">tc_grps</A></b>)
-is set with <b><A HREF="#tc">ref_t</A></b> [K],
+(<b><A HREF="#tc">tc-grps</A></b>)
+is set with <b><A HREF="#tc">ref-t</A></b> [K],
the inverse friction constant for each group is set with
-<b><A HREF="#tc">tau_t</A></b> [ps].
+<b><A HREF="#tc">tau-t</A></b> [ps].
The parameter <b><A HREF="#tc">tcoupl</A></b> is ignored.
-The random generator is initialized with <b><A HREF="#ld">ld_seed</A></b>.
-When used as a thermostat, an appropriate value for <b>tau_t</b> is 2 ps,
+The random generator is initialized with <b><A HREF="#ld">ld-seed</A></b>.
+When used as a thermostat, an appropriate value for <b>tau-t</b> is 2 ps,
since this results in a friction that is lower than the internal friction
of water, while it is high enough to remove excess heat
(unless <b>cut-off</b> or <b>reaction-field</b> electrostatics is used).
NOTE: temperature deviations decay twice as fast as with
-a Berendsen thermostat with the same <b>tau_t</b>.</dd>
+a Berendsen thermostat with the same <b>tau-t</b>.</dd>
<dt><b>sd1</b></dt>
<dd> An efficient leap-frog stochastic dynamics integrator.
This integrator is equivalent to <b>sd</b>, except that it requires
-only one Gaussian random number and one constraint step.
-This integrator is less accurate.
-For water the relative error in the temperature with this integrator
-is 0.5 <b>delta_t</b>/<b>tau_t</b>, but for other systems it can be higher.
-Use with care and check the temperature.</dd>
+only one Gaussian random number and one constraint step and is therefore
+significantly faster. Without constraints the accuracy is the same as <b>sd</b>.
+With constraints the accuracy is significantly reduced, so then <b>sd</b>
+will often be preferred.</dd>
<dt><b>bd</b></dt>
<dd>An Euler integrator for Brownian or position Langevin dynamics, the
velocity is the force divided by a friction coefficient
-(<b><A HREF="#ld">bd_fric</A></b> [amu ps<sup>-1</sup>])
-plus random thermal noise (<b><A HREF="#tc">ref_t</A></b>).
-When <b><A HREF="#ld">bd_fric</A></b><tt>=0</tt>, the friction coefficient for each
-particle is calculated as mass/<b><A HREF="#tc">tau_t</A></b>, as for the
+(<b><A HREF="#ld">bd-fric</A></b> [amu ps<sup>-1</sup>])
+plus random thermal noise (<b><A HREF="#tc">ref-t</A></b>).
+When <b><A HREF="#ld">bd-fric</A></b><tt>=0</tt>, the friction coefficient for each
+particle is calculated as mass/<b><A HREF="#tc">tau-t</A></b>, as for the
integrator <tt>sd</tt>.
-The random generator is initialized with <b><A HREF="#ld">ld_seed</A></b>.</dd>
+The random generator is initialized with <b><A HREF="#ld">ld-seed</A></b>.</dd>
<dt><b>steep</b></dt>
<dd>A <!--Idx-->steepest descent<!--EIdx--> algorithm for energy
which is only allowed for single-atom molecules).
Since neighborlist construction is expensive, one can perform several
extra insertions with the same list almost for free.
-The random seed is set with <b><A HREF="#ld">ld_seed</A></b>.
+The random seed is set with <b><A HREF="#ld">ld-seed</A></b>.
The temperature for the Boltzmann weighting is set with
-<b><A HREF="#tc">ref_t</A></b>, this should match the temperature
+<b><A HREF="#tc">ref-t</A></b>, this should match the temperature
of the simulation of the original trajectory.
Dispersion correction is implemented correctly for tpi.
All relevant quantities are written to the file specified with
<tt>sd</tt> and <tt>bd</tt>)</dd>
<dt><b>nsteps: (0)</b></dt>
<dd>maximum number of steps to integrate or minimize, -1 is no maximum</dd>
-<dt><b>init_step: (0)</b></dt>
+<dt><b>init-step: (0)</b></dt>
<dd>The starting step.
-The time at an step i in a run is calculated as: t = <tt>tinit</tt> + <tt>dt</tt>*(<tt>init_step</tt> + i).
-The free-energy lambda is calculated as: lambda = <tt>init_lambda</tt> + <tt>delta_lambda</tt>*(<tt>init_step</tt> + i).
+The time at an step i in a run is calculated as: t = <tt>tinit</tt> + <tt>dt</tt>*(<tt>init-step</tt> + i).
+The free-energy lambda is calculated as: lambda = <tt>init-lambda</tt> + <tt>delta-lambda</tt>*(<tt>init-step</tt> + i).
Also non-equilibrium MD parameters can depend on the step number.
Thus for exact restarts or redoing part of a run it might be necessary to
-set <tt>init_step</tt> to the step number of the restart frame.
+set <tt>init-step</tt> to the step number of the restart frame.
<tt>tpbconv</tt> does this automatically.
</dd>
-<dt><b>comm_mode:</b></dt>
+<dt><b>comm-mode:</b></dt>
<dd><dl compact>
<dt><b>Linear</b></dt>
<dd>Remove center of mass translation</dd>
<dt><b>Angular</b></dt>
<dd>Remove center of mass translation and rotation around the center of mass
</dd>
- <dt><b>No</b></dt>
+ <dt><b>None</b></dt>
<dd>No restriction on the center of mass motion
</dl></dd>
<dt><b>nstcomm: (10) [steps]</b></dt>
<dd>frequency for center of mass motion removal</dd>
-<dt><b>comm_grps:</b></dt>
+<dt><b>comm-grps:</b></dt>
<dd>group(s) for center of mass motion removal, default is the whole system</dd>
</dl>
<h3><!--Idx-->Langevin dynamics<!--EIdx--></h3>
<dl>
-<dt><b>bd_fric: (0) [amu ps<sup>-1</sup>]</b></dt>
+<dt><b>bd-fric: (0) [amu ps<sup>-1</sup>]</b></dt>
<dd>Brownian dynamics friction coefficient.
-When <b>bd_fric</b><tt>=0</tt>, the friction coefficient for each
-particle is calculated as mass/<b><A HREF="#tc">tau_t</A></b>.</dd>
-<dt><b>ld_seed: (1993) [integer]</b></dt>
+When <b>bd-fric</b><tt>=0</tt>, the friction coefficient for each
+particle is calculated as mass/<b><A HREF="#tc">tau-t</A></b>.</dd>
+<dt><b>ld-seed: (1993) [integer]</b></dt>
<dd>used to initialize random generator for thermal noise
for stochastic and Brownian dynamics.
-When <b>ld_seed</b> is set to -1, the seed is calculated from the process ID.
+When <b>ld-seed</b> is set to -1, the seed is calculated from the process ID.
When running BD or SD on multiple processors, each processor uses a seed equal
-to <b>ld_seed</b> plus the processor number.</dd>
+to <b>ld-seed</b> plus the processor number.</dd>
</dl>
<A NAME="em"><br>
<hr>
<h3>Output control</h3>
<dl>
-<dt><b>nstxout: (100) [steps]</b></dt>
+<dt><b>nstxout: (0) [steps]</b></dt>
<dd>frequency to write coordinates to output
<!--Idx-->trajectory file<!--EIdx-->, the last coordinates are always written</dd>
-<dt><b>nstvout: (100) [steps]</b></dt>
+<dt><b>nstvout: (0) [steps]</b></dt>
<dd>frequency to write velocities to output trajectory,
the last velocities are always written</dd>
<dt><b>nstfout: (0) [steps]</b></dt>
<dd>frequency to write forces to output trajectory.</dd>
-<dt><b>nstlog: (100) [steps]</b></dt>
+<dt><b>nstlog: (1000) [steps]</b></dt>
<dd>frequency to write energies to <!--Idx-->log file<!--EIdx-->,
the last energies are always written</dd>
<dt><b>nstcalcenergy: (-1)</b></dt>
energy averages and fluctuations also when <b>nstenergy</b><tt>>1</tt></dd>
<dt><b>nstxtcout: (0) [steps]</b></dt>
<dd>frequency to write coordinates to xtc trajectory</dd>
-<dt><b>xtc_precision: (1000) [real]</b></dt>
+<dt><b>xtc-precision: (1000) [real]</b></dt>
<dd>precision to write to xtc trajectory</dd>
-<dt><b>xtc_grps:</b></dt>
+<dt><b>xtc-grps:</b></dt>
<dd>group(s) to write to xtc trajectory, default the whole system is written
(if <b>nstxtcout</b> > 0)</dd>
<dt><b>energygrps:</b></dt>
</dd>
</dl></dd>
-<dt><b>ns_type:</b></dt>
+<dt><b>ns-type:</b></dt>
<dd><dl compact>
<dt><b>grid</b></dt>
<dd>Make a grid in the box and only check atoms in neighboring grid
<dd>Use no periodic boundary conditions, ignore the box.
To simulate without cut-offs, set all cut-offs to 0 and <b>nstlist</b><tt>=0</tt>.
For best performance without cut-offs, use <b>nstlist</b><tt>=0</tt>,
-<b>ns_type</b><tt>=simple</tt>
+<b>ns-type</b><tt>=simple</tt>
and particle decomposition instead of domain decomposition.</dd>
<dt><b>xy</b></dt>
<dd>Use periodic boundary conditions in x and y directions only.
-This works only with <b>ns_type</b><tt>=grid</tt> and can be used
+This works only with <b>ns-type</b><tt>=grid</tt> and can be used
in combination with <b><a href="#walls">walls</a></b>.
Without walls or with only one wall the system size is infinite
in the z direction. Therefore pressure coupling or Ewald summation
These disadvantages do not apply when two walls are used.</dd>
</dl></dd>
-<dt><b>periodic_molecules:</b></dt>
+<dt><b>periodic-molecules:</b></dt>
<dd><dl compact>
<dt><b>no</b></dt>
<dd>molecules are finite, fast molecular PBC can be used</dd>
and molecules are not made whole in the output</dd>
</dl></dd>
-<dt><b>rlist: (1) [nm]</b></dt>
-<dd>cut-off distance for the short-range neighbor list</dd>
+<dt><b>rlist: (-1) [nm]</b></dt>
+<dd>cut-off distance for the short-range neighbor list, should be ≥ 0</dd>
<dt><b>rlistlong: (-1) [nm]</b></dt>
<dd>Cut-off distance for the long-range neighbor list.
Use e.g. <b>rlist</b><tt>=0.9</tt>, <b>rcoulomb</b><tt>=0.9</tt>. The highest magnitude of
wave vectors used in reciprocal space is controlled by <b>fourierspacing</b>.
The relative accuracy of direct/reciprocal space
-is controlled by <b>ewald_rtol</b>.
+is controlled by <b>ewald-rtol</b>.
<br>
NOTE: Ewald scales as O(N<sup>3/2</sup>)
and is thus extremely slow for large systems. It is included mainly for
reference - in most cases PME will perform much better.</dd>
<dt><b><!--Idx-->PME<!--EIdx--></b></dt>
-<dd>Fast Particle-Mesh Ewald electrostatics. Direct space is similar
+<dd>Fast smooth Particle-Mesh Ewald (SPME) electrostatics. Direct space is similar
to the Ewald sum, while the reciprocal part is performed with
FFTs. Grid dimensions are controlled with <b>fourierspacing</b> and the
-interpolation order with <b>pme_order</b>. With a grid spacing of 0.1
+interpolation order with <b>pme-order</b>. With a grid spacing of 0.1
nm and cubic interpolation the electrostatic forces have an accuracy
of 2-3*10<sup>-4</sup>. Since the error from the vdw-cutoff is larger than this you
might try 0.15 nm. When running in parallel the interpolation
parallelizes better than the FFT, so try decreasing grid dimensions
while increasing interpolation.</dd>
-<dt><b><!--Idx-->PPPM<!--EIdx--></b></dt>
-<dd>Particle-Particle Particle-Mesh algorithm for long range
-electrostatic interactions.
-Use for example <b>rlist</b><tt>=0.9</tt>, <b>rcoulomb</b><tt>=0.9</TT>.
-The grid dimensions are controlled by <b>fourierspacing</b>.
-Reasonable grid spacing for PPPM is 0.05-0.1 nm.
-See <tt>Shift</tt> for the details of the particle-particle potential.
-<br>
-NOTE: PPPM is not functional in the current version, but we plan to implement
-PPPM through a small modification of the PME code.</dd>
+<dt><b><!--Idx-->P3M-AD<!--EIdx--></b></dt>
+<dd>Particle-Particle Particle-Mesh algorithm with analytical derivative
+for for long range electrostatic interactions. The method and code
+is identical to SPME, except that the influence function is optimized
+for the grid. This gives a slight increase in accuracy.</dd>
<dt><b>Reaction-Field electrostatics<!--QuietIdx-->reaction-field electrostatics<!--EQuietIdx--></b></dt>
<dd>Reaction field with Coulomb cut-off <b>rcoulomb</b>,
where <b>rcoulomb</b> ≥ <b>rlist</b>.
-The dielectric constant beyond the cut-off is <b>epsilon_rf</b>.
-The dielectric constant can be set to infinity by setting <b>epsilon_rf</b><tt>=0</tt>.</dd>
+The dielectric constant beyond the cut-off is <b>epsilon-rf</b>.
+The dielectric constant can be set to infinity by setting <b>epsilon-rf</b><tt>=0</tt>.</dd>
<dt><b>Generalized-Reaction-Field</b></dt>
<dd>Generalized reaction field with Coulomb cut-off <b>rcoulomb</b>,
where <b>rcoulomb</b> ≥ <b>rlist</b>.
-The dielectric constant beyond the cut-off is <b>epsilon_rf</b>.
+The dielectric constant beyond the cut-off is <b>epsilon-rf</b>.
The ionic strength is computed from the number of charged
(i.e. with non zero charge) <!--Idx-->charge group<!--EIdx-->s.
The temperature for the GRF potential is set with
-<b><A HREF="#tc">ref_t</A></b> [K].</dd>
+<b><A HREF="#tc">ref-t</A></b> [K].</dd>
<dt><b>Reaction-Field-zero</b></dt>
<dd>In GROMACS normal reaction-field electrostatics leads to bad
energy conservation. <b>Reaction-Field-zero</b> solves this
by making the potential zero beyond the cut-off. It can only
-be used with an infinite dielectric constant (<b>epsilon_rf=0</b>),
+be used with an infinite dielectric constant (<b>epsilon-rf=0</b>),
because only for that value the force vanishes at the cut-off.
<b>rlist</b> should be 0.1 to 0.3 nm larger than <b>rcoulomb</b>
to accommodate for the size of charge groups and diffusion
</dl></dd>
<A NAME="el2">
-<dt><b>rcoulomb_switch: (0) [nm]</b></dt>
+<dt><b>rcoulomb-switch: (0) [nm]</b></dt>
<dd>where to start switching the Coulomb potential</dd>
-<dt><b>rcoulomb: (1) [nm]</b></dt>
-<dd>distance for the Coulomb <!--Idx-->cut-off<!--EIdx--></dd>
+<dt><b>rcoulomb: (-1) [nm]</b></dt>
+<dd>distance for the Coulomb <!--Idx-->cut-off<!--EIdx-->, should be ≥ 0</dd>
-<dt><b>epsilon_r: (1)</b></dt>
+<dt><b>epsilon-r: (1)</b></dt>
<dd>The relative <!--Idx-->dielectric constant<!--EIdx-->.
A value of 0 means infinity.</dd>
-<dt><b>epsilon_rf: (1)</b></dt>
+<dt><b>epsilon-rf: (0)</b></dt>
<dd>The relative dielectric constant of the reaction field.
This is only used with reaction-field electrostatics.
A value of 0 means infinity.</dd>
where <b>rvdw</b> <tt>≥</tt> <b>rlist</b>.</dd>
<dt><b>Shift</b></dt>
<dd>The LJ (not Buckingham) potential is decreased over the whole
-range and the forces decay smoothly to zero between <b>rvdw_switch</b>
+range and the forces decay smoothly to zero between <b>rvdw-switch</b>
and <b>rvdw</b>. The neighbor search cut-off <b>rlist</b> should be
0.1 to 0.3 nm larger than <b>rvdw</b> to accommodate for the size of
charge groups and diffusion between neighbor list
<dt><b>Switch</b></dt>
<dd>The LJ (not Buckingham)
-potential is normal out to <b>rvdw_switch</b>, after which it is switched
+potential is normal out to <b>rvdw-switch</b>, after which it is switched
off to reach zero at <b>rvdw</b>. Both the potential and force functions
are continuously smooth, but be aware that all switch functions will give rise
to a bulge (increase) in the force (since we are switching the potential).
for <tt>f</tt> and <tt>-f'</tt> are ignored.</dd>
</dl></dd>
-<dt><b>rvdw_switch: (0) [nm]</b></dt>
+<dt><b>rvdw-switch: (0) [nm]</b></dt>
<dd>where to start switching the LJ potential</dd>
-<dt><b>rvdw: (1) [nm]</b></dt>
-<dd>distance for the LJ or Buckingham <!--Idx-->cut-off<!--EIdx--></dd>
+<dt><b>rvdw: (-1) [nm]</b></dt>
+<dd>distance for the LJ or Buckingham <!--Idx-->cut-off<!--EIdx-->, should be ≥ 0</dd>
<dt><b>DispCorr:</b></dt>
<dd><dl compact></dd>
which are always tabulated irrespective of the use of
tables for the non-bonded interactions. </dd>
-<dt><b>energygrp_table:</b></dt>
+<dt><b>energygrp-table:</b></dt>
<dd>When user tables are used for electrostatics and/or VdW,
here one can give pairs of energy groups for which seperate
user tables should be used.
The two energy groups will be appended to the table file name,
in order of their definition in <b>energygrps</b>, seperated by underscores.
For example, if <tt>energygrps = Na Cl Sol</tt>
-and <tt>energygrp_table = Na Na Na Cl</tt>, <tt>mdrun</tt> will read
+and <tt>energygrp-table = Na Na Na Cl</tt>, <tt>mdrun</tt> will read
<tt>table_Na_Na.xvg</tt> and <tt>table_Na_Cl.xvg</tt> in addition
to the normal <tt>table.xvg</tt> which will be used for all other
energy group pairs.
<dt><b>fourierspacing: (0.12) [nm]</b></dt>
<dd>For ordinary Ewald, the ratio of the box dimensions and the spacing
determines a lower bound for the number of wave vectors to use in each
-(signed) direction. For PME and PPPM, that ratio determines a lower bound
+(signed) direction. For PME and P3M, that ratio determines a lower bound
for the number of Fourier-space grid points that will be used along that
axis. In all cases, the number for each direction can be overridden by
entering a non-zero value for <b>fourier_n[xyz]</b>.
when the Coulomb cut-off and the PME grid spacing are scaled
by the same factor.</dd>
-<dt><b>fourier_nx (0) ; fourier_ny (0) ; fourier_nz: (0)</b></dt>
+<dt><b>fourier-nx (0) ; fourier-ny (0) ; fourier-nz: (0)</b></dt>
<dd>Highest magnitude of wave vectors in reciprocal space when using Ewald.</dd>
-<dd>Grid size when using PPPM or PME. These values override
+<dd>Grid size when using PME or P3M. These values override
<b>fourierspacing</b> per direction. The best choice is powers of
2, 3, 5 and 7. Avoid large primes.</dd>
-<dt><b>pme_order (4)</b></dt>
+<dt><b>pme-order (4)</b></dt>
<dd>Interpolation order for PME. 4 equals cubic interpolation. You might try
6/8/10 when running in parallel and simultaneously decrease grid dimension.</dd>
-<dt><b>ewald_rtol (1e-5)</b></dt>
+<dt><b>ewald-rtol (1e-5)</b></dt>
<dd>The relative strength of the Ewald-shifted direct potential at
-<b>rcoulomb</b> is given by <b>ewald_rtol</b>.
+<b>rcoulomb</b> is given by <b>ewald-rtol</b>.
Decreasing this will give a more accurate direct sum,
but then you need more wave vectors for the reciprocal sum.</dd>
-<dt><b>ewald_geometry: (3d)</b></dt>
+<dt><b>ewald-geometry: (3d)</b></dt>
<dd><dl compact>
<dt><b>3d</b></dt>
<dd>The Ewald sum is performed in all three dimensions.</dd>
and use this option.</dd>
</dl></dd>
-<dt><b>epsilon_surface: (0)</b></dt>
+<dt><b>epsilon-surface: (0)</b></dt>
<dd>This controls the dipole correction to the Ewald summation in 3D. The
default value of zero means it is turned off. Turn it on by setting it to the value
of the relative permittivity of the imaginary surface around your infinite system. Be
This value does not affect the slab 3DC variant of the long range corrections.</dd>
-<dt><b>optimize_fft:</b></dt>
+<dt><b>optimize-fft:</b></dt>
<dd><dl compact>
<dt><b>no</b></dt>
<dd>Don't calculate the optimal FFT plan for the grid at startup.</dd>
<dd>No temperature coupling.</dd>
<dt><b>berendsen</b></dt>
<dd>Temperature coupling with a Berendsen-thermostat to a bath with
-temperature <b>ref_t</b> [K], with time constant <b>tau_t</b> [ps].
+temperature <b>ref-t</b> [K], with time constant <b>tau-t</b> [ps].
Several groups can be coupled separately, these are specified in the
-<b>tc_grps</b> field separated by spaces.</dd>
+<b>tc-grps</b> field separated by spaces.</dd>
<dt><b>nose-hoover</b></dt>
<dd>Temperature coupling using a Nose-Hoover extended
ensemble. The reference temperature and coupling groups are selected
-as above, but in this case <b>tau_t</b> [ps] controls the period
+as above, but in this case <b>tau-t</b> [ps] controls the period
of the temperature fluctuations at equilibrium, which is slightly
different from a relaxation time.
For NVT simulations the conserved energy quantity is written
<dd>Temperature coupling using velocity rescaling with a stochastic term
(JCP 126, 014101).
This thermostat is similar to Berendsen coupling, with the same scaling
-using <b>tau_t</b>, but the stochastic term ensures that a proper
+using <b>tau-t</b>, but the stochastic term ensures that a proper
canonical ensemble is generated. The random seed is set with
-<b><A HREF="#ld">ld_seed</A></b>.
-This thermostat works correctly even for <b>tau_t</b><tt>=0</tt>.
+<b><A HREF="#ld">ld-seed</A></b>.
+This thermostat works correctly even for <b>tau-t</b><tt>=0</tt>.
For NVT simulations the conserved energy quantity is written
to the energy and log file.</dd>
</dl>
</dd>
<dt><b>nh-chain-length (10)</b></dt>
<dd>the number of chained Nose-Hoover thermostats for velocity Verlet integrators, the leap-frog <b>md</b> integrator only supports 1. Data for the NH chain variables is not printed to the .edr, but can be using the <tt>GMX_NOSEHOOVER_CHAINS</tt> environment variable</dd>
-<dt><b>tc_grps:</b></dt>
+<dt><b>tc-grps:</b></dt>
<dd>groups to couple separately to temperature bath</dd>
-<dt><b>tau_t: [ps]</b></dt>
-<dd>time constant for coupling (one for each group in <b>tc_grps</b>),
+<dt><b>tau-t: [ps]</b></dt>
+<dd>time constant for coupling (one for each group in <b>tc-grps</b>),
-1 means no temperature coupling</dd>
-<dt><b>ref_t: [K]</b></dt>
-<dd>reference temperature for coupling (one for each group in <b>tc_grps</b>)</dd>
+<dt><b>ref-t: [K]</b></dt>
+<dd>reference temperature for coupling (one for each group in <b>tc-grps</b>)</dd>
</dl>
<A NAME="pc"><br>
<dd>No pressure coupling. This means a fixed box size.</dd>
<dt><b>berendsen</b></dt>
<dd>Exponential relaxation pressure coupling with time constant
-<b>tau_p</b> [ps]. The box is scaled every timestep. It has been
+<b>tau-p</b> [ps]. The box is scaled every timestep. It has been
argued that this does not yield a correct thermodynamic ensemble,
but it is the most efficient way to scale a box at the beginning
of a run.</dd>
<dd>Extended-ensemble pressure coupling where the box vectors are
subject to an equation of motion. The equation of motion for the atoms
is coupled to this. No instantaneous scaling takes place. As for
-Nose-Hoover temperature coupling the time constant <b>tau_p</b> [ps]
+Nose-Hoover temperature coupling the time constant <b>tau-p</b> [ps]
is the period of pressure fluctuations at equilibrium. This is
probably a better method when you want to apply pressure scaling
during data collection, but beware that you can get very large
<dt><b>MTTK</b></dt>
<dd>Martyna-Tuckerman-Tobias-Klein implementation, only useable with <b>md-vv</b>
or <b>md-vv-avek</b>, very similar to Parrinello-Rahman.
-As for Nose-Hoover temperature coupling the time constant <b>tau_p</b>
+As for Nose-Hoover temperature coupling the time constant <b>tau-p</b>
[ps] is the period of pressure fluctuations at equilibrium. This is
probably a better method when you want to apply pressure scaling
during data collection, but beware that you can get very large
<dt><b>pcoupltype:</b></dt>
<dd><dl compact>
<dt><b>isotropic</b></dt>
-<dd>Isotropic pressure coupling with time constant <b>tau_p</b> [ps].
+<dd>Isotropic pressure coupling with time constant <b>tau-p</b> [ps].
The compressibility and reference pressure are set with
-<b>compressibility</b> [bar<sup>-1</sup>] and <b>ref_p</b> [bar], one
+<b>compressibility</b> [bar<sup>-1</sup>] and <b>ref-p</b> [bar], one
value is needed.</dd>
<dt><b>semiisotropic</b></dt>
<dd>Pressure coupling which is isotropic in the <tt>x</tt> and <tt>y</tt> direction,
<dd>Surface tension coupling for surfaces parallel to the xy-plane.
Uses normal pressure coupling for the <tt>z</tt>-direction, while the surface tension
is coupled to the <tt>x/y</tt> dimensions of the box.
-The first <b>ref_p</b> value is the reference surface tension times
+The first <b>ref-p</b> value is the reference surface tension times
the number of surfaces [bar nm],
the second value is the reference <tt>z</tt>-pressure [bar].
The two <b>compressibility</b> [bar<sup>-1</sup>] values are the compressibility
For velocity Verlet integrators <b>nstpcouple</b> is set to 1.</dd>
</dd>
-<dt><b>tau_p: (1) [ps]</b></dt>
+<dt><b>tau-p: (1) [ps]</b></dt>
<dd>time constant for coupling</dd>
<dt><b>compressibility: [bar<sup>-1</sup>]</b></dt>
<dd>compressibility (NOTE: this is now really in bar<sup>-1</sup>)
For water at 1 atm and 300 K the compressibility is 4.5e-5 [bar<sup>-1</sup>].</dd>
-<dt><b>ref_p: [bar]</b></dt>
+<dt><b>ref-p: [bar]</b></dt>
<dd>reference pressure for coupling</dd>
-<dt><b>refcoord_scaling:</b></dt>
+<dt><b>refcoord-scaling:</b></dt>
<dd><dl compact>
<dt><b>no</b></dt>
<dd>The reference coordinates for position restraints are not modified.
</dd>
</dl>
-<dt><b>annealing_npoints:</b></dt>
+<dt><b>annealing-npoints:</b></dt>
<dd>A list with the number of annealing reference/control points used for
each temperature group. Use 0 for groups that are not annealed. The number of entries should equal the number of temperature groups.</dd>
-<dt><b>annealing_time:</b></dt>
-<dd>List of times at the annealing reference/control points for each group. If you are using periodic annealing, the times will be used modulo the last value, i.e. if the values are 0, 5, 10, and 15, the coupling will restart at the 0ps value after 15ps, 30ps, 45ps, etc. The number of entries should equal the sum of the numbers given in <tt>annealing_npoints</tt>.</dd>
+<dt><b>annealing-time:</b></dt>
+<dd>List of times at the annealing reference/control points for each group. If you are using periodic annealing, the times will be used modulo the last value, i.e. if the values are 0, 5, 10, and 15, the coupling will restart at the 0ps value after 15ps, 30ps, 45ps, etc. The number of entries should equal the sum of the numbers given in <tt>annealing-npoints</tt>.</dd>
-<dt><b>annealing_temp:</b></dt>
-<dd>List of temperatures at the annealing reference/control points for each group. The number of entries should equal the sum of the numbers given in <tt>annealing_npoints</tt>.</dd>
+<dt><b>annealing-temp:</b></dt>
+<dd>List of temperatures at the annealing reference/control points for each group. The number of entries should equal the sum of the numbers given in <tt>annealing-npoints</tt>.</dd>
<br>
-Confused? OK, let's use an example. Assume you have two temperature groups, set the group selections to <tt>annealing = single periodic</tt>, the number of points of each group to <tt>annealing_npoints = 3 4</tt>, the times to <tt>annealing_time = 0 3 6 0 2 4 6</tt> and finally temperatures to <tt>annealing_temp = 298 280 270 298 320 320 298</tt>.
+Confused? OK, let's use an example. Assume you have two temperature groups, set the group selections to <tt>annealing = single periodic</tt>, the number of points of each group to <tt>annealing-npoints = 3 4</tt>, the times to <tt>annealing-time = 0 3 6 0 2 4 6</tt> and finally temperatures to <tt>annealing-temp = 298 280 270 298 320 320 298</tt>.
The first group will be coupled to 298K at 0ps, but the reference temperature will drop linearly to reach 280K at 3ps, and then linearly between 280K and 270K from 3ps to 6ps. After this is stays constant, at 270K. The second group is coupled to 298K at 0ps, it increases linearly to 320K at 2ps, where it stays constant until 4ps. Between 4ps and 6ps it decreases to 298K, and then it starts over with the same pattern again, i.e. rising linearly from 298K to 320K between 6ps and 8ps. Check the summary printed by <tt>grompp</tt> if you are unsure!
</dl>
<h3>Velocity generation</h3>
<dl>
-<dt><b>gen_vel:</b></dt>
+<dt><b>gen-vel:</b></dt>
<dd><dl compact>
<dt><b>no</b></dt>
<dd> Do not generate velocities. The velocities are set to zero
when there are no velocities in the input structure file.</dd>
<dt><b>yes</b></dt>
<dd>Generate velocities in <tt>grompp</tt> according to a Maxwell distribution at
-temperature <b>gen_temp</b> [K], with random seed <b>gen_seed</b>.
+temperature <b>gen-temp</b> [K], with random seed <b>gen-seed</b>.
This is only meaningful with integrator <b><A HREF="#run">md</A></b>.</dd>
</dl></dd>
-<dt><b>gen_temp: (300) [K]</b></dt>
+<dt><b>gen-temp: (300) [K]</b></dt>
<dd>temperature for Maxwell distribution</dd>
-<dt><b>gen_seed: (173529) [integer]</b></dt>
+<dt><b>gen-seed: (173529) [integer]</b></dt>
<dd>used to initialize random generator for random velocities,
-when <b>gen_seed</b> is set to -1, the seed is calculated from
+when <b>gen-seed</b> is set to -1, the seed is calculated from
the process ID number.
</dl>
<dd>Convert all bonds and angles to bond-constraints.</dd>
</dl>
-<dt><b>constraint_algorithm:</b></dt>
+<dt><b>constraint-algorithm:</b></dt>
<dd><dl compact>
<dt><b><!--Idx-->LINCS<!--EIdx--></b></dt>
<dd>LINear Constraint Solver.
With domain decomposition the parallel version P-LINCS is used.
The accuracy in set with
-<b>lincs_order</b>, which sets the number of matrices in the expansion
+<b>lincs-order</b>, which sets the number of matrices in the expansion
for the matrix inversion.
After the matrix inversion correction the algorithm does
an iterative correction to compensate for lengthening due to rotation.
The number of such iterations can be controlled with
-<b>lincs_iter</b>. The root mean square relative constraint deviation
+<b>lincs-iter</b>. The root mean square relative constraint deviation
is printed to the log file every <b>nstlog</b> steps.
-If a bond rotates more than <b>lincs_warnangle</b> [degrees] in one step,
+If a bond rotates more than <b>lincs-warnangle</b> [degrees] in one step,
a warning will be printed both to the log file and to <TT>stderr</TT>.
LINCS should not be used with coupled angle constraints.
</dd>
<dt><b><!--Idx-->SHAKE<!--EIdx--></b></dt>
<dd>SHAKE is slightly slower and less stable than LINCS, but does work with
angle constraints.
-The relative tolerance is set with <b>shake_tol</b>, 0.0001 is a good value
+The relative tolerance is set with <b>shake-tol</b>, 0.0001 is a good value
for ``normal'' MD. SHAKE does not support constraints between atoms
on different nodes, thus it can not be used with domain decompositon
when inter charge-group constraints are present.
</dd>
</dl></dd>
<dt><b>continuation:</b></dt>
-<dd>This option was formerly known as <tt>unconstrained_start</tt>.</dd>
+<dd>This option was formerly known as <tt>unconstrained-start</tt>.</dd>
<dd><dl compact>
<dt><b>no</b></dt>
<dd>apply constraints to the start configuration and reset shells</dd>
</dl></dd>
<A NAME="bond2">
-<dt><b>shake_tol: (0.0001)</b></dt>
+<dt><b>shake-tol: (0.0001)</b></dt>
<dd>relative tolerance for SHAKE</dd>
-<dt><b>lincs_order: (4)</b></dt>
+<dt><b>lincs-order: (4)</b></dt>
<dd>Highest order in the expansion of the constraint coupling matrix.
When constraints form triangles, an additional expansion of the same
order is applied on top of the normal expansion only for the couplings
needed for large time-steps with virtual sites or BD.
For accurate energy minimization an order of 8 or more might be required.
With domain decomposition, the cell size is limited by the distance
-spanned by <b>lincs_order</b>+1 constraints. When one wants to scale
-further than this limit, one can decrease <b>lincs_order</b> and increase
-<b>lincs_iter</b>, since the accuracy does not deteriorate
-when (1+<b>lincs_iter</b>)*<b>lincs_order</b> remains constant.</dd>
-<dt><b>lincs_iter: (1)</b></dt>
+spanned by <b>lincs-order</b>+1 constraints. When one wants to scale
+further than this limit, one can decrease <b>lincs-order</b> and increase
+<b>lincs-iter</b>, since the accuracy does not deteriorate
+when (1+<b>lincs-iter</b>)*<b>lincs-order</b> remains constant.</dd>
+<dt><b>lincs-iter: (1)</b></dt>
<dd>Number of iterations to correct for rotational lengthening in LINCS.
For normal runs a single step is sufficient, but for NVE
runs where you want to conserve energy accurately or for accurate
energy minimization you might want to increase it to 2.
-<dt><b>lincs_warnangle: </b>(30) [degrees]</dt>
+<dt><b>lincs-warnangle: </b>(30) [degrees]</dt>
<dd>maximum angle that a bond can rotate before LINCS will complain</dd>
<dt><b>morse:</b></dt>
<hr>
<h3>Energy group <!--Idx-->exclusions<!--EIdx--></h3>
<dl>
-<dt><b>energygrp_excl: </b></dt>
+<dt><b>energygrp-excl: </b></dt>
<dd>Pairs of energy groups for which all non-bonded interactions are
excluded. An example: if you have two energy groups <tt>Protein</tt>
and <tt>SOL</tt>, specifying
<br>
-<tt>energygrp_excl = Protein Protein SOL SOL</tt>
+<tt>energygrp-excl = Protein Protein SOL SOL</tt>
<br>
would give only the non-bonded interactions between the protein and the
solvent. This is especially useful for speeding up energy calculations with
<dl>
<dt><b>nwall: 0</b></dt>
<dd>When set to <b>1</b> there is a wall at <tt>z=0</tt>, when set to <b>2</b>
-there is also a wall at <tt>z=z_box</tt>. Walls can only be used with <b>pbc=xy</b>.
+there is also a wall at <tt>z=z-box</tt>. Walls can only be used with <b>pbc=xy</b>.
When set to <b>2</b> pressure coupling and Ewald summation can be used
(it is usually best to use semiisotropic pressure coupling with
the <tt>x/y</tt> compressibility set to 0, as otherwise the surface area will change).
-Walls interact wit the rest of the system through an optional <tt>wall_atomtype</tt>.
+Walls interact wit the rest of the system through an optional <tt>wall-atomtype</tt>.
Energy groups <tt>wall0</tt> and <tt>wall1</tt> (for <b>nwall=2</b>) are
added automatically to monitor the interaction of energy groups
with each wall.
The <A HREF="#run">center of mass motion removal</A> will be turned
off in the <tt>z</tt>-direction.</dd>
-<dt><b>wall_atomtype:</b></dt>
+<dt><b>wall-atomtype:</b></dt>
<dd>the atom type name in the force field for each wall.
By (for example) defining a special wall atom type in the topology with its
own combination rules, this allows for independent tuning of the interaction
of each atomtype with the walls.</dd>
-<dt><b>wall_type:</b></dt>
+<dt><b>wall-type:</b></dt>
<dd><dl compact>
<dt><b>9-3</b></dt>
<dd>LJ integrated over the volume behind the wall: 9-3 potential</dd>
<dt><b>12-6</b></dt>
<dd>direct LJ potential with the z distance from the wall</dd>
<dt><b>table</b></dt><dd>user defined potentials indexed with the z distance from the wall, the tables are read analogously to
-the <b><A HREF="#table">energygrp_table</A></b> option,
+the <b><A HREF="#table">energygrp-table</A></b> option,
where the first name is for a ``normal'' energy group and the second name
is <tt>wall0</tt> or <tt>wall1</tt>,
only the dispersion and repulsion columns are used</dd>
</dl></dd>
-<dt><b>wall_r_linpot: -1 (nm)</b></dt>
+<dt><b>wall-r-linpot: -1 (nm)</b></dt>
<dd>Below this distance from the wall the potential is continued
linearly and thus the force is constant. Setting this option to
a postive value is especially useful for equilibration when some atoms
are beyond a wall.
-When the value is ≤0 (<0 for <b>wall_type=table</b>),
+When the value is ≤0 (<0 for <b>wall-type=table</b>),
a fatal error is generated when atoms are beyond a wall.
</dd>
-<dt><b>wall_density: [nm<sup>-3</sup>/nm<sup>-2</sup>]</b></dt>
+<dt><b>wall-density: [nm<sup>-3</sup>/nm<sup>-2</sup>]</b></dt>
<dd>the number density of the atoms for each wall for wall types
<b>9-3</b> and <b>10-4</b>
-<dt><b>wall_ewald_zfac: 3</b></dt>
+<dt><b>wall-ewald-zfac: 3</b></dt>
<dd>The scaling factor for the third box vector for Ewald summation only,
the minimum is 2.
Ewald summation can only be used with <b>nwall=2</b>, where one
-should use <b><A HREF="#ewald">ewald_geometry</A><tt>=3dc</tt></b>.
+should use <b><A HREF="#ewald">ewald-geometry</A><tt>=3dc</tt></b>.
The empty layer in the box serves to decrease the unphysical Coulomb
interaction between periodic images.
</dl>
between the reference group and one or more groups.
The setup is identical to the option <b>umbrella</b>, except for the fact
that a rigid constraint is applied instead of a harmonic potential.</dd>
-<dt><b>constant_force</b></dt>
+<dt><b>constant-force</b></dt>
<dd>Center of mass pulling using a linear potential and therefore
a constant force. For this option there is no reference position
-and therefore the parameters <b>pull_init</b> and <b>pull_rate</b>
+and therefore the parameters <b>pull-init</b> and <b>pull-rate</b>
are not used.</dd>
</dl></dd>
-<dt><b>pull_geometry:</b></dt>
+<dt><b>pull-geometry:</b></dt>
<dd><dl compact>
<dt><b>distance</b></dt>
<dd>Pull along the vector connecting the two groups.
-Components can be selected with <b>pull_dim</b>.</dd>
+Components can be selected with <b>pull-dim</b>.</dd>
<dt><b>direction</b></dt>
-<dd>Pull in the direction of <b>pull_vec</b>.</dd>
-<dt><b>direction_periodic</b></dt>
+<dd>Pull in the direction of <b>pull-vec</b>.</dd>
+<dt><b>direction-periodic</b></dt>
<dd>As <b>direction</b>, but allows the distance to be larger than
half the box size. With this geometry the box should not be dynamic
(e.g. no pressure scaling) in the pull dimensions and the pull force
<dt><b>cylinder</b></dt>
<dd>Designed for pulling with respect to a layer where the reference COM
is given by a local cylindrical part of the reference group.
-The pulling is in the direction of <b>pull_vec</b>.
+The pulling is in the direction of <b>pull-vec</b>.
From the reference group a cylinder is selected around the axis going
-through the pull group with direction <b>pull_vec</b> using two radii.
-The radius <b>pull_r1</b> gives the radius within which all
-the relative weights are one, between <b>pull_r1</b> and
-<b>pull_r0</b> the weights are switched to zero. Mass weighting is also used.
+through the pull group with direction <b>pull-vec</b> using two radii.
+The radius <b>pull-r1</b> gives the radius within which all
+the relative weights are one, between <b>pull-r1</b> and
+<b>pull-r0</b> the weights are switched to zero. Mass weighting is also used.
Note that the radii should be smaller than half the box size.
For tilted cylinders they should be even smaller than half the box size
since the distance of an atom in the reference group
from the COM of the pull group has both a radial and an axial component.
<dt><b>position</b></dt>
<dd>Pull to the position of the reference group plus
-<b>pull_init</b> + time*<b>pull_rate</b>*<b>pull_vec</b>.</dd>
+<b>pull-init</b> + time*<b>pull-rate</b>*<b>pull-vec</b>.</dd>
</dl></dd>
-<dt><b>pull_dim: (Y Y Y)</b></dt>
+<dt><b>pull-dim: (Y Y Y)</b></dt>
<dd>the distance components to be used with geometry <b>distance</b>
and <b>position</b>, and also sets which components are printed
to the output files</dd>
-<dt><b>pull_r1: (1) [nm]</b></dt>
+<dt><b>pull-r1: (1) [nm]</b></dt>
<dd>the inner radius of the cylinder for geometry <b>cylinder</b></dd>
-<dt><b>pull_r0: (1) [nm]</b></dt>
+<dt><b>pull-r0: (1) [nm]</b></dt>
<dd>the outer radius of the cylinder for geometry <b>cylinder</b></dd>
-<dt><b>pull_constr_tol: (1e-6)</b></dt>
+<dt><b>pull-constr-tol: (1e-6)</b></dt>
<dd>the relative constraint tolerance for constraint pulling</dd>
-<dt><b>pull_start:</b></dt>
+<dt><b>pull-start:</b></dt>
<dd><dl compact>
<dt><b>no</b></dt>
-<dd>do not modify <b>pull_init</b>
+<dd>do not modify <b>pull-init</b>
<dt><b>yes</b></dt>
-<dd>add the COM distance of the starting conformation to <b>pull_init</b></dd>
+<dd>add the COM distance of the starting conformation to <b>pull-init</b></dd>
</dl>
-<dt><b>pull_nstxout: (10)</b></dt>
+<dt><b>pull-nstxout: (10)</b></dt>
<dd>frequency for writing out the COMs of all the pull group</dd>
-<dt><b>pull_nstfout: (1)</b></dt>
+<dt><b>pull-nstfout: (1)</b></dt>
<dd>frequency for writing out the force of all the pulled group</dd>
-<dt><b>pull_ngroups: (1)</b></dt>
+<dt><b>pull-ngroups: (1)</b></dt>
<dd>The number of pull groups, not including the reference group.
If there is only one group, there is no difference in treatment
of the reference and pulled group (except with the cylinder geometry).
and the first group (ending on 1) are given,
further groups work analogously, but with the number 1 replaced
by the group number.</dd>
-<dt><b>pull_group0: </b></dt>
+<dt><b>pull-group0: </b></dt>
<dd>The name of the reference group. When this is empty an absolute reference
of (0,0,0) is used. With an absolute reference the system is no longer
translation invariant and one should think about what to do with
the <A HREF="#run">center of mass motion</A>.</dd>
-<dt><b>pull_weights0: </b></dt>
-<dd>see <b>pull_weights1</b></dd>
-<dt><b>pull_pbcatom0: (0)</b></dt>
-<dd>see <b>pull_pbcatom1</b></dd>
-<dt><b>pull_group1: </b></dt>
+<dt><b>pull-weights0: </b></dt>
+<dd>see <b>pull-weights1</b></dd>
+<dt><b>pull-pbcatom0: (0)</b></dt>
+<dd>see <b>pull-pbcatom1</b></dd>
+<dt><b>pull-group1: </b></dt>
<dd>The name of the pull group.</dd>
-<dt><b>pull_weights1: </b></dt>
+<dt><b>pull-weights1: </b></dt>
<dd>Optional relative weights which are multiplied with the masses of the atoms
to give the total weight for the COM. The number should be 0, meaning all 1,
or the number of atoms in the pull group.</dd>
-<dt><b>pull_pbcatom1: (0)</b></dt>
+<dt><b>pull-pbcatom1: (0)</b></dt>
<dd>The reference atom for the treatment of periodic boundary conditions
inside the group
(this has no effect on the treatment of the pbc between groups).
This option is only important when the diameter of the pull group
is larger than half the shortest box vector.
For determining the COM, all atoms in the group are put at their periodic image
-which is closest to <b>pull_pbcatom1</b>.
+which is closest to <b>pull-pbcatom1</b>.
A value of 0 means that the middle atom (number wise) is used.
This parameter is not used with geometry <b>cylinder</b>.
A value of -1 turns on cosine weighting, which is useful for a group
of molecules in a periodic system, e.g. a water slab (see Engin et al.
J. Chem. Phys. B 2010).</dd>
-<dt><b>pull_vec1: (0.0 0.0 0.0)</b></dt>
+<dt><b>pull-vec1: (0.0 0.0 0.0)</b></dt>
<dd>The pull direction. <tt>grompp</tt> normalizes the vector.</dd>
-<dt><b>pull_init1: (0.0) / (0.0 0.0 0.0) [nm]</b></dt>
+<dt><b>pull-init1: (0.0) / (0.0 0.0 0.0) [nm]</b></dt>
<dd>The reference distance at t=0. This is a single value,
except for geometry <b>position</b> which uses a vector.</dd>
-<dt><b>pull_rate1: (0) [nm/ps]</b></dt>
+<dt><b>pull-rate1: (0) [nm/ps]</b></dt>
<dd>The rate of change of the reference position.</dd>
-<dt><b>pull_k1: (0) [kJ mol<sup>-1</sup> nm<sup>-2</sup>] / [kJ mol<sup>-1</sup> nm<sup>-1</sup>]</b></dt>
+<dt><b>pull-k1: (0) [kJ mol<sup>-1</sup> nm<sup>-2</sup>] / [kJ mol<sup>-1</sup> nm<sup>-1</sup>]</b></dt>
<dd>The force constant. For umbrella pulling this is the harmonic force
constant in [kJ mol<sup>-1</sup> nm<sup>-2</sup>]. For constant force pulling
this is the force constant of the linear potential, and thus minus (!)
the constant force in [kJ mol<sup>-1</sup> nm<sup>-1</sup>].</dd>
-<dt><b>pull_kB1: (pull_k1) [kJ mol<sup>-1</sup> nm<sup>-2</sup>] / [kJ mol<sup>-1</sup> nm<sup>-1</sup>]</b></dt>
-<dd>As <b>pull_k1</b>, but for state B. This is only used when
-<A HREF="#free"><b>free_energy</b></A> is turned on.
-The force constant is then (1 - lambda)*<b>pull_k1</b> + lambda*<b>pull_kB1</b>.
+<dt><b>pull-kB1: (pull-k1) [kJ mol<sup>-1</sup> nm<sup>-2</sup>] / [kJ mol<sup>-1</sup> nm<sup>-1</sup>]</b></dt>
+<dd>As <b>pull-k1</b>, but for state B. This is only used when
+<A HREF="#free"><b>free-energy</b></A> is turned on.
+The force constant is then (1 - lambda)*<b>pull-k1</b> + lambda*<b>pull-kB1</b>.
</dl>
<A NAME="nmr"><br>
of systems within each ensemble (usually equal to the <tt>mdrun -multi</tt> value).</dd>
</dd>
</dl></dd>
-<dt><b>disre_weighting:</b></dt>
+<dt><b>disre-weighting:</b></dt>
<dd><dl compact>
<dt><b>equal</b> (default)</dt>
<dd>divide the restraint force equally over all atom pairs in the restraint</dd>
<dt><b>conservative</b></dt>
<dd>the forces are the derivative of the restraint potential,
this results in an r<sup>-7</sup> weighting of the atom pairs.
-The forces are conservative when <tt>disre_tau</tt> is zero.</dd>
+The forces are conservative when <tt>disre-tau</tt> is zero.</dd>
</dl></dd>
-<dt><b>disre_mixed:</b></dt>
+<dt><b>disre-mixed:</b></dt>
<dd><dl compact>
<dt><b>no</b></dt>
<dd>the violation used in the calculation of the restraint force is the
square root of the product of the time-averaged violation and the instantaneous violation</dd>
</dl></dd>
-<dt><b>disre_fc: (1000) [kJ mol<sup>-1</sup> nm<sup>-2</sup>]</b></dt>
+<dt><b>disre-fc: (1000) [kJ mol<sup>-1</sup> nm<sup>-2</sup>]</b></dt>
<dd>force constant for distance restraints, which is multiplied by a
(possibly) different factor for each restraint given in the <tt>fac</tt>
column of the interaction in the topology file.</dd>
-<dt><b>disre_tau: (0) [ps]</b></dt>
+<dt><b>disre-tau: (0) [ps]</b></dt>
<dd>time constant for distance restraints running average. A value of zero turns off time averaging.</dd>
<dt><b>nstdisreout: (100) [steps]</b></dt>
<dd>use orientation restraints, ensemble averaging can be performed
with <tt>mdrun -multi</tt></dd>
</dl>
-<dt><b>orire_fc: (0) [kJ mol]</b></dt>
+<dt><b>orire-fc: (0) [kJ mol]</b></dt>
<dd>force constant for orientation restraints, which is multiplied by a
(possibly) different weight factor for each restraint, can be set to zero to
obtain the orientations from a free simulation</dd>
-<dt><b>orire_tau: (0) [ps]</b></dt>
+<dt><b>orire-tau: (0) [ps]</b></dt>
<dd>time constant for orientation restraints running average. A value of zero turns off time averaging.</dd>
-<dt><b>orire_fitgrp: </b></dt>
+<dt><b>orire-fitgrp: </b></dt>
<dd>fit group for orientation restraining. This group of atoms is used
to determine the rotation <b>R</b> of the system with respect to the
reference orientation. The reference orientation is the starting
<h3>Free energy calculations<!--QuietIdx-->free energy calculations<!--EQuietIdx--></h3>
<dl>
-<dt><b>free_energy:</b></dt>
+<dt><b>free-energy:</b></dt>
<dd><dl compact>
<dt><b>no</b></dt>
<dd>Only use topology A.</dd>
<dt><b>yes</b></dt>
<dd>Interpolate between topology A (lambda=0) to topology B (lambda=1)
-and write the derivative of the Hamiltonian with respect to lambda (as specified with <b>dhdl_derivatives</b>), or the Hamiltonian differences with respect to other lambda values (as specified with <b>foreign_lambda</b>) to
+and write the derivative of the Hamiltonian with respect to lambda (as specified with <b>dhdl-derivatives</b>), or the Hamiltonian differences with respect to other lambda values (as specified with <b>foreign-lambda</b>) to
the energy file and/or to <tt>dhdl.xvg</tt>, where they can be processed by, for example <tt>g_bar</tt>.
The potentials, bond-lengths and angles are interpolated linearly as
-described in the manual. When <b>sc_alpha</b> is larger than zero, soft-core
+described in the manual. When <b>sc-alpha</b> is larger than zero, soft-core
potentials are used for the LJ and Coulomb interactions.</dd>
</dl></dd>
-<dt><b>init_lambda: (0)</b></dt>
-<dd>starting value for lambda</dd>
-<dt><b>delta_lambda: (0)</b></dt>
+<dt><b>init-lambda: (0)</b></dt>
+<dd>starting value for lambda (float). Generally, this should only be used with slow growth. In other cases, <b>init-lambda-state</b> should be specified instead.</dd>
+<dt><b>init-lambda-state: (0)</b></dt>
+<dd>starting value for the lambda state (integer). Specified which columm of the lambda vector should be used.</dd>
+<dt><b>delta-lambda: (0)</b></dt>
<dd>increment per time step for lambda</dd>
-<dt><b>foreign_lambda: ()</b></dt>
+<dt><b>coul-lambdas: ()</b></dt>
+<dd>Zero, one or more lambda values for which Delta H values will
+be determined and written to dhdl.xvg every <b>nstdhdl</b> steps.
+Only the electrostatic interactions are controlled with this component of the lambda vector.</dd>
+<dt><b>vdw-lambdas: ()</b></dt>
+<dd>Zero, one or more lambda values for which Delta H values will
+be determined and written to dhdl.xvg every <b>nstdhdl</b> steps.
+Only the van der Waals interactions are controlled with this component of the lambda vector.</dd>
+<dt><b>bonded-lambdas: ()</b></dt>
+<dd>Zero, one or more lambda values for which Delta H values will
+be determined and written to dhdl.xvg every <b>nstdhdl</b> steps.
+Only the bonded interactions are controlled with this component of the lambda vector.</dd>
+<dt><b>restraint-lambdas: ()</b></dt>
+<dd>Zero, one or more lambda values for which Delta H values will
+be determined and written to dhdl.xvg every <b>nstdhdl</b> steps.
+Only the restraint interactions are controlled with this component of the lambda vector.</dd>
+<dt><b>mass-lambdas: ()</b></dt>
+<dd>Zero, one or more lambda values for which Delta H values will
+be determined and written to dhdl.xvg every <b>nstdhdl</b> steps.
+Only the particle masses are controlled with this component of the lambda vector.</dd>
+<dt><b>temperature-lambdas: ()</b></dt>
+<dd>Zero, one or more lambda values for which Delta H values will
+be determined and written to dhdl.xvg every <b>nstdhdl</b> steps.
+Only the temperatures controlled with this component of the lambda vector.
+Note that these lambdas should not be used for replica exchange, only for simulated tempering.</dd>
+<dt><b>fep-lambdas: ()</b></dt>
<dd>Zero, one or more lambda values for which Delta H values will
be determined and written to dhdl.xvg every <b>nstdhdl</b> steps.
Free energy differences between different lambda values can then
-be determined with <tt>g_bar</tt>.</dd>
-<dt><b>dhdl_derivatives: (yes)</b></dt>
+be determined with <tt>g_bar</tt>. <b>fep-lambdas</b> is different from the other -lambdas keywords because
+all components of the lambda vector that are not specified will use <b>fep-lambdas</b>.</dd>
+<dt><b>dhdl-derivatives: (yes)</b></dt>
<dd>If yes (the default), the derivatives of the Hamiltonian with respect to lambda at each <b>nstdhdl</b> step are written out. These values are needed for interpolation of linear energy differences with <tt>g_bar</tt> (although the same can also be achieved with the right <b>foreign lambda</b> setting, that may not be as flexible), or with thermodynamic integration</dd>
-<dt><b>sc_alpha: (0)</b></dt>
+<dt><b>sc-alpha: (0)</b></dt>
<dd>the soft-core parameter, a value of 0 results in linear interpolation of
the LJ and Coulomb interactions</dd>
-<dt><b>sc_power: (0)</b></dt>
+<dt><b>sc-r-power: (6)</b></dt>
+<dd>the power of the radial term in the soft-core equation. Possible values are 6 and 48. 6 is more standard, and is the default. When 48 is used, then sc-alpha should generally be much lower (between 0.001 and 0.003).</dd>
+<dt><b>sc-coul: (no)</b></dt>
+<dd>Whether to apply the soft core free energy interations to the Columbic interaction. Default is no, as it is generally
+more efficient to turn of the Coulomic interactions linearly before turning off electrostatic interactions.</dd>
+<dt><b>sc-power: (0)</b></dt>
<dd>the power for lambda in the soft-core function,
only the values 1 and 2 are supported</dd>
-<dt><b>sc_sigma: (0.3) [nm]</b></dt>
+<dt><b>sc-sigma: (0.3) [nm]</b></dt>
<dd>the soft-core sigma for particles which have a C6 or C12 parameter equal
-to zero or a sigma smaller than <b>sc_sigma</b></dd>
+to zero or a sigma smaller than <b>sc-sigma</b></dd>
<dt><b>couple-moltype:</b></dt>
<dd>Here one can supply a molecule type (as defined in the topology)
for calculating solvation or coupling free energies.
There is a special option <b>system</b> that couples all molecule types
in the system. This can be useful for equilibrating a system
starting from (nearly) random coordinates.
-<b>free_energy</b> has to be turned on.
+<b>free-energy</b> has to be turned on.
The Van der Waals interactions and/or charges in this molecule type can be
turned on or off between lambda=0 and lambda=1, depending on the settings
of <b>couple-lambda0</b> and <b>couple-lambda1</b>. If you want to decouple
</dl>
<dt><b>nstdhdl: (10)</b></dt>
<dd>the frequency for writing dH/dlambda and possibly Delta H to dhdl.xvg,
-0 means no ouput, should be a multiple of <b>nstcalcenergy</b></dd>
-<dt><b>separate_dhdl_file: (yes)</b></dt>
+0 means no ouput, should be a multiple of <b>nstcalcenergy</b>and <b>nstfep</b></dd>
+<dt><b>nstfep: (10)</b></dt>
+<dd>the frequency at which energies at other values of lambda are calculated. If not specified, set to be the same as <b>nstdhdl</b>. Should be a multiple of <b>nstcalcenergy</b>. If replica exchange is chosen, then -replex must also be a multiple of <b>nstfep</b></dd>
+<dt><b>separate-dhdl-file: (yes)</b></dt>
<dd><dl compact>
<dt><b>yes</b></dt>
-<dd>the free energy values that are calculated (as specified with the <b>foreign-lambda</b> and <b>dhdl_derivatives</b> settings) are written out to a separate file, with the default name <tt>dhdl.xvg</tt>. This file can be used directly with <tt>g_bar</tt>.</dd>
+<dd>the free energy values that are calculated (as specified with the <b>foreign-lambda</b> and <b>dhdl-derivatives</b> settings) are written out to a separate file, with the default name <tt>dhdl.xvg</tt>. This file can be used directly with <tt>g_bar</tt>.</dd>
<dt><b>no</b></dt>
<dd>The free energy values are written out to the energy output file (<tt>ener.edr</tt>, in accumulated blocks at every <b>nstenergy</b> steps), where they can be extracted with <tt>g_energy</tt> or used directly with <tt>g_bar</tt>.</dd>
</dl>
-<dt><b>dh_hist_size: (0)</b></dt>
-<dd>If nonzero, specifies the size of the histogram into which the Delta H values (specified with <b>foreign_lambda</b>) and the derivative dH/dl values are binned, and written to ener.edr. This can be used to save disk space while calculating free energy differences. One histogram gets written for each <b>foreign lambda</b> and two for the dH/dl, at every <b>nstenergy</b> step. Be aware that incorrect histogram settings (too small size or too wide bins) can introduce errors. Do not use histograms unless you're certain you need it.</dd>
-<dt><b>dh_hist_spacing (0.1)</b></dt>
-<dd>Specifies the bin width of the histograms, in energy units. Used in conjunction with <b>dh_hist_size</b>. This size limits the accuracy with which free energies can be calculated. Do not use histograms unless you're certain you need it.</dd>
+<dt><b>dh-hist-size: (0)</b></dt>
+<dd>If nonzero, specifies the size of the histogram into which the Delta H values (specified with <b>foreign-lambda</b>) and the derivative dH/dl values are binned, and written to ener.edr. This can be used to save disk space while calculating free energy differences. One histogram gets written for each <b>foreign lambda</b> and two for the dH/dl, at every <b>nstenergy</b> step. Be aware that incorrect histogram settings (too small size or too wide bins) can introduce errors. Do not use histograms unless you're certain you need it.</dd>
+<dt><b>dh-hist-spacing (0.1)</b></dt>
+<dd>Specifies the bin width of the histograms, in energy units. Used in conjunction with <b>dh-hist-size</b>. This size limits the accuracy with which free energies can be calculated. Do not use histograms unless you're certain you need it.</dd>
+</dl>
+<A NAME="expanded"><br>
+<hr>
+<h3><!--Idx-->Expanded Ensemble calculations<!--EIdx--></h3>
+
+<dl>
+<dt><b>nstexpanded</b></dt> <dd>The frequency to peform expanded ensemble
+simulations. Must be a multiple of <b>nstfep</b>.</dd>
+<dt><b>lmc-stats:</b></dt>
+<dd><dl compact>
+<dt><b>no</b></dt>
+<dd>No Monte Carlo in state space</dd>
+<dt><b>metropolis-transition</b></dt>
+<dd> Uses the Metropolis weights to update the expanded ensemble weight of the state.
+Min{1,exp(-(beta_new u_new - beta_old u_old)}</dd>
+<dt><b>barker-transition</b></dt>
+<dd> Uses the Barker transition critera to update the expanded ensemble weight of the state.</dd>
+<dt><b>wang-landau</b></dt>
+<dd>Uses the Wang-Landau algorithm (in state space) to update the expanded ensemble weights.</dd>
+<dt><b>min-variance</b></dt>
+<dd>Uses the minimum variance updating method of Escobedo et al to update the expanded ensemble weights. Weights
+will not be the free energies, but will rather emphasize states that need more sampling to give even uncertainty.
+</dl>
+<dt><b>lmc-mc-move:</b></dt>
+<dd><dl compact>
+<dt><b>no</b></dt>
+<dd>No Monte Carlo in state space is performed.</dd>
+<dt><b>metropolis-transition</b></dt>
+<dd> Randomly chooses a new state up or down, then uses the Metropolis critera to decide whether to accept or reject:
+Min{1,exp(-(beta_new u_new - beta_old u_old)}</dd>
+<dt><b>barker-transition</b></dt>
+<dd> Randomly chooses a new state up or down, then uses the Barker transition critera to decide whether to accept or reject: exp(-beta_new u_new)/[exp(-beta_new u_new)+exp(-beta_old u_old)] </dd>
+<dt><b>gibbs</b></dt>
+<dd> Uses the conditional weights of the state given the coordinate (exp(-beta_i u_i) / sum_k exp(beta_i u_i) to
+decide which state to move to.</dd>
+<dt><b>metropolized-gibbs</b></dt>
+<dd>
+<dd> Uses the conditional weights of the state given the coordinate (exp(-beta_i u_i) / sum_k exp(beta_i u_i) to
+decide which state to move to, EXCLUDING the current state, then uses a rejection step to ensure detailed
+balance. Always more efficient that Gibbs, though marginally so in many situations.</dd>
</dl>
+<dt><b>lmc-seed:</b></dt>
+<dd> random seed to use for Monte Carlo moves in state space. If not specified, <b>ld-seed</b> is used instead. </dd>
+<dt><b>mc-temperature:</b></dt>
+<dd> Temperature used for acceptance/rejection for Monte Carlo moves. If not specified, the temperature of the
+simulation specified in the first group of <b>ref_t</b> is used.</dd>
+
+<dt><b>wl-scale: (0.8)</b></dt>
+<dt><b>wl-ratio: (0.8)</b></dt>
+<dt><b>init-wl-delta: (1.0) </b></dt>
+<dt><b>wl-oneovert: (no) </b></dt>
+<dt><b>lmc-repeats: (1)</b></dt>
+<dt><b>lmc-gibbsdelta: (-1) </b></dt>
+<dt><b>lmc-forced-nstart: (0) </b></dt>
+<dt><b>nst-transition-matrix: (-1)</b></dt>
+<dd>Frequency of outputting the expanded ensemble transition matrix. A negative number means it will only be printed at the end of the simulation.<dd>
+<dt><b>symmetrized-transition-matrix: (no) </b></dt>
+<dd>Whether to symmetrize the empirical transition matrix</dd>
+<dt><b>mininum-var-min</b></dt>
+<dt><b>weight-c-range</b></dt>
+
+<dt><b>simulated-tempering: (no)</b></dt>
+<dt><b>simulated-tempering-scaling: ()</b></dt>
+<dt><b>sim-temp-low: (300):</b></dt>
+<dd>Low temperature for simulated tempering</dd>
+<dt><b>sim-temp-high: (300):</b></dt>
+<dd>High temperature for simulated tempering</dd>
+</dl>
<A NAME="neq"><br>
<hr>
<h3>Non-equilibrium MD<!--QuietIdx-->non-equilibrium MD<!--EQuietIdx--></h3>
<dl>
-<dt><b>acc_grps: </b></dt>
+<dt><b>acc-grps: </b></dt>
<dd>groups for constant acceleration (e.g.: <tt>Protein Sol</tt>)
all atoms in groups Protein and Sol will experience constant acceleration
as specified in the <b>accelerate</b> line</dd>
<dt><b>accelerate: (0) [nm ps<sup>-2</sup>]</b></dt>
-<dd>acceleration for <b>acc_grps</b>; x, y and z for each group
+<dd>acceleration for <b>acc-grps</b>; x, y and z for each group
(e.g. <tt>0.1 0.0 0.0 -0.1 0.0 0.0</tt> means that first group has constant
acceleration of 0.1 nm ps<sup>-2</sup> in X direction, second group the
opposite).</dd>
(e.g. <tt>Y Y N N N N</tt> means that particles in the first group
can move only in Z direction. The particles in the second group can
move in any direction).</dd>
-<dt><b>cos_acceleration: (0) [nm ps<sup>-2</sup>]</b></dt>
+<dt><b>cos-acceleration: (0) [nm ps<sup>-2</sup>]</b></dt>
<dd>the amplitude of the acceleration profile for calculating the
<!--Idx-->viscosity<!--EIdx-->.
The acceleration is in the X-direction and the magnitude is
-<b>cos_acceleration</b> cos(2 pi z/boxheight).
+<b>cos-acceleration</b> cos(2 pi z/boxheight).
Two terms are added to the energy file:
the amplitude of the velocity profile and 1/viscosity.</dd>
<dt><b><!--Idx-->deform<!--EIdx-->: (0 0 0 0 0 0) [nm ps<sup>-1</sup>]</b></dt>
<h3>Electric fields<!--QuietIdx-->electric field<!--EQuietIdx--></h3>
<dl>
-<dt><b>E_x ; E_y ; E_z:</b></dt>
+<dt><b>E-x ; E-y ; E-z:</b></dt>
<dd>If you want to use an electric field in a direction, enter 3 numbers
-after the appropriate <b>E_*</b>, the first number: the number of cosines,
+after the appropriate <b>E-*</b>, the first number: the number of cosines,
only 1 is implemented (with frequency 0) so enter 1,
the second number: the strength of the electric field in
<b>V nm<sup>-1</sup></b>,
the third number: the phase of the cosine, you can enter any number here
since a cosine of frequency zero has no phase.</dd>
-<dt><b>E_xt; E_yt; E_zt: </b></dt>
+<dt><b>E-xt; E-yt; E-zt: </b></dt>
<dd>not implemented yet</dd>
</dl>
<br>
<h3>Implicit solvent</h3>
<dl>
-<dt><b>implicit_solvent:</b></dt>
+<dt><b>implicit-solvent:</b></dt>
<dd><dl compact="compact">
<dt><b>no</b></dt>
<dd>No implicit solvent</dd>
<dt><b>GBSA</b></dt>
<dd>Do a simulation with implicit solvent using the Generalized Born formalism.
Three different methods for calculating the Born radii are available, Still, HCT and
-OBC. These are specified with the <b>gb_algorithm</b> field. The non-polar solvation
-is specified with the <b>sa_algorithm</b> field.</dd>
+OBC. These are specified with the <b>gb-algorithm</b> field. The non-polar solvation
+is specified with the <b>sa-algorithm</b> field.</dd>
</dl>
-<dt><b>gb_algorithm:</b></dt>
+<dt><b>gb-algorithm:</b></dt>
<dd><dl compact="compact">
<dt><b>Still</b></dt>
<dd>Use the Still method to calculate the Born radii</dd>
<dt><b>rgbradii: (1.0) [nm]</b></dt>
<dd>Cut-off for the calculation of the Born radii. Currently must be equal to rlist</dd>
-<dt><b>gb_epsilon_solvent: (80)</b></dt>
+<dt><b>gb-epsilon-solvent: (80)</b></dt>
<dd>Dielectric constant for the implicit solvent</dd>
-<dt><b>gb_saltconc: (0) [M]</b></dt>
+<dt><b>gb-saltconc: (0) [M]</b></dt>
<dd>Salt concentration for implicit solvent models, currently not used</dd>
-<dt><b>gb_obc_alpha (1); gb_obc_beta (0.8); gb_obc_gamma (4.85);</b></dt>
+<dt><b>gb-obc-alpha (1); gb-obc-beta (0.8); gb-obc-gamma (4.85);</b></dt>
<dd>Scale factors for the OBC model. Default values are OBC(II).
Values for OBC(I) are 0.8, 0 and 2.91 respectively</dd>
-<dt><b>gb_dielectric_offset: (0.009) [nm]</b></dt>
+<dt><b>gb-dielectric-offset: (0.009) [nm]</b></dt>
<dd>Distance for the di-electric offset when calculating the Born radii. This is
the offset between the center of each atom the center of the polarization energy
for the corresponding atom</dd>
-<dt><b>sa_algorithm</b></dt>
+<dt><b>sa-algorithm</b></dt>
<dd><dl compact="compact">
<dt><b>Ace-approximation</b></dt>
<dd>Use an Ace-type approximation (default)</dd>
calculated</dd>
</dl>
-<dt><b>sa_surface_tension: [kJ mol<sup>-1</sup> nm<sup>-2</sup>]</b></dt>
+<dt><b>sa-surface-tension: [kJ mol<sup>-1</sup> nm<sup>-2</sup>]</b></dt>
<dd>Default value for surface tension with SA algorithms. The default value is -1;
Note that if this default value is not changed
it will be overridden by <tt>grompp</tt> using values that are specific for the choice
of radii algorithm (0.0049 kcal/mol/Angstrom<sup>2</sup> for Still, 0.0054 kcal/mol/Angstrom<sup>2</sup>
for HCT/OBC)
-Setting it to 0 will while using an sa_algorithm other than None means
+Setting it to 0 will while using an sa-algorithm other than None means
no non-polar calculations are done.
</dd>
</dl>
+<A NAME="adress"><br>
+<hr>
+<h3>Adaptive Resolution Simulation</h3>
+
+<dl>
+<dt><b>adress: (no)</b></dt>
+<dd>Decide whether the AdResS feature is turned on.</dd>
+<dt><b>adress-type: (Off)</b></dt>
+<dd><dl compact>
+<dt><b>Off</b></dt>
+<dd>Do an AdResS simulation with weight equal 1, which is equivalent to an explicit (normal) MD simulation. The difference to disabled AdResS is that the AdResS variables are still read-in and hence are defined.</dd>
+<dt><b>Constant</b></dt>
+<dd>Do an AdResS simulation with a constant weight, <b>adress-const-wf</b> defines the value of the weight</dd>
+<dt><b>XSplit</b></dt>
+<dd>Do an AdResS simulation with simulation box split in x-direction, so basically the weight is only a function of the x coordinate and all distances are measured using the x coordinate only.</dd>
+<dt><b>Sphere</b></dt>
+<dd>Do an AdResS simulation with spherical explicit zone.</dd>
+</dl></dd>
+<dt><b>adress-const-wf: (1)</b></dt>
+<dd>Provides the weight for a constant weight simulation (<b>adress-type</b>=Constant)</dd>
+<dt><b>adress-ex-width: (0)</b></dt>
+<dd>Width of the explicit zone, measured from <b>adress-reference-coords</b>.</dd>
+<dt><b>adress-hy-width: (0)</b></dt>
+<dd>Width of the hybrid zone.</dd>
+<dt><b>adress-reference-coords: (0,0,0)</b></dt>
+<dd>Position of the center of the explicit zone. Periodic boundary conditions apply for measuring the distance from it.</dd>
+<dt><b>adress-cg-grp-names</b></dt>
+<dd>The names of the coarse-grained energy groups. All other energy groups are considered explicit and their interactions will be automatically excluded with the coarse-grained groups.</dd>
+<dt><b>adress-site: (COM)</b>The mapping point from which the weight is calculated.</dt>
+<dd><dl compact>
+<dt><b>COM</b></dt>
+<dd>The weight is calculated from the center of mass of each charge group.</dd>
+<dt><b>COG</b></dt>
+<dd>The weight is calculated from the center of geometry of each charge group.</dd>
+<dt><b>Atom</b></dt>
+<dd>The weight is calculated from the position of 1st atom of each charge group.</dd>
+<dt><b>AtomPerAtom</b></dt>
+<dd>The weight is calculated from the position of each individual atom.</dd>
+</dl></dd>
+<dt><b>adress-interface-correction: (Off)</b></dt>
+<dd><dl compact>
+<dt><b>Off</b></dt>
+<dd>Do not a apply any interface correction.</dd>
+<dt><b>thermoforce</b></dt>
+<dd>Apply thermodynamic force interface correction. The table can be specified using the <tt>-tabletf</tt> option of <tt>mdrun</tt>. The table should contain the potential and force (acting on molecules) as function of the distance from <b>adress-reference-coords</b>.</dd>
+</dl></dd>
+<dt><b>adress-tf-grp-names</b></dt>
+<dd>The names of the energy groups to which the <b>thermoforce</b> is applied if enabled in <b>adress-interface-correction</b>. If no group is given the default table is applied.</dd>
+<dt><b>adress-ex-forcecap: (0)</b></dt>
+<dd>Cap the force in the hybrid region, useful for big molecules. 0 disables force capping.</dd>
+</dl>
+
<A NAME="user"><br>
<hr>
<h3>User defined thingies</h3>
<dl>
-<dt><b>user1_grps; user2_grps: </b></dt>
+<dt><b>user1-grps; user2-grps: </b></dt>
<dt><b>userint1 (0); userint2 (0); userint3 (0); userint4 (0)</b></dt>
<dt><b>userreal1 (0); userreal2 (0); userreal3 (0); userreal4 (0)</b></dt>
<dd>These you can use if you modify code. You can pass integers and
<P>
<multicol cols=4>
-<A HREF="#neq">acc_grps</A><br>
+<A HREF="#neq">acc-grps</A><br>
<A HREF="#neq">accelerate</A><br>
<A HREF="#sa">annealing</A><br>
-<A HREF="#sa">annealing_npoints</A><br>
-<A HREF="#sa">annealing_time</A><br>
-<A HREF="#sa">annealing_temp</A><br>
-<A HREF="#ld">bd_fric</A><br>
+<A HREF="#sa">annealing-npoints</A><br>
+<A HREF="#sa">annealing-time</A><br>
+<A HREF="#sa">annealing-temp</A><br>
+<A HREF="#ld">bd-fric</A><br>
<A HREF="#vdw">bDispCorr</A><br>
-<A HREF="#run">comm_mode</A><br>
-<A HREF="#run">comm_grps</A><br>
+<A HREF="#run">comm-mode</A><br>
+<A HREF="#run">comm-grps</A><br>
<A HREF="#pc">compressibility</A><br>
-<A HREF="#bond">constraint_algorithm</A><br>
+<A HREF="#bond">constraint-algorithm</A><br>
<A HREF="#bond">constraints</A><br>
-<A HREF="#neq">cos_acceleration</A><br>
+<A HREF="#neq">cos-acceleration</A><br>
<A HREF="#el">coulombtype</A><br>
<A HREF="#free">couple-intramol</A><br>
<A HREF="#free">couple-lambda0</A><br>
<A HREF="#free">couple-moltype</A><br>
<A HREF="#pp">define</A><br>
<A HREF="#neq">deform</A><br>
-<A HREF="#free">delta_lambda</A><br>
+<A HREF="#free">delta-lambda</A><br>
<A HREF="#nmr">disre</A><br>
-<A HREF="#nmr">disre_weighting</A><br>
-<A HREF="#nmr">disre_mixed</A><br>
-<A HREF="#nmr">disre_fc</A><br>
-<A HREF="#nmr">disre_tau</A><br>
+<A HREF="#nmr">disre-weighting</A><br>
+<A HREF="#nmr">disre-mixed</A><br>
+<A HREF="#nmr">disre-fc</A><br>
+<A HREF="#nmr">disre-tau</A><br>
<A HREF="#run">dt</A><br>
<A HREF="#em">emstep</A><br>
<A HREF="#em">emtol</A><br>
-<A HREF="#egexcl">energygrp_excl</A><br>
-<A HREF="#table">energygrp_table</A><br>
+<A HREF="#egexcl">energygrp-excl</A><br>
+<A HREF="#table">energygrp-table</A><br>
<A HREF="#out">energygrps</A><br>
-<A HREF="#el2">epsilon_r</A><br>
-<A HREF="#el2">epsilon_rf</A><br>
-<A HREF="#ewald">ewald_rtol</A><br>
-<A HREF="#ewald">ewald_geometry</A><br>
-<A HREF="#ewald">epsilon_surface</A><br>
-<A HREF="#ef">E_x</A><br>
-<A HREF="#ef">E_xt</A><br>
-<A HREF="#ef">E_y</A><br>
-<A HREF="#ef">E_yt</A><br>
-<A HREF="#ef">E_z</A><br>
-<A HREF="#ef">E_zt </A><br>
+<A HREF="#el2">epsilon-r</A><br>
+<A HREF="#el2">epsilon-rf</A><br>
+<A HREF="#ewald">ewald-rtol</A><br>
+<A HREF="#ewald">ewald-geometry</A><br>
+<A HREF="#ewald">epsilon-surface</A><br>
+<A HREF="#ef">E-x</A><br>
+<A HREF="#ef">E-xt</A><br>
+<A HREF="#ef">E-y</A><br>
+<A HREF="#ef">E-yt</A><br>
+<A HREF="#ef">E-z</A><br>
+<A HREF="#ef">E-zt </A><br>
<A HREF="#xmdrun">fcstep</A><br>
-<A HREF="#ewald">fourier_nx</A><br>
-<A HREF="#ewald">fourier_ny</A><br>
-<A HREF="#ewald">fourier_nz</A><br>
+<A HREF="#ewald">fourier-nx</A><br>
+<A HREF="#ewald">fourier-ny</A><br>
+<A HREF="#ewald">fourier-nz</A><br>
<A HREF="#ewald">fourierspacing</A><br>
-<A HREF="#free">free_energy</A><br>
+<A HREF="#free">free-energy</A><br>
<A HREF="#neq">freezedim </A><br>
<A HREF="#neq">freezegrps</A><br>
-<A HREF="#vel">gen_seed</A><br>
-<A HREF="#vel">gen_temp</A><br>
-<A HREF="#vel">gen_vel</A><br>
+<A HREF="#vel">gen-seed</A><br>
+<A HREF="#vel">gen-temp</A><br>
+<A HREF="#vel">gen-vel</A><br>
<A HREF="#pp">include</A><br>
-<A HREF="#free">init_lambda</A><br>
-<A HREF="#run">init_step</A><br>
+<A HREF="#free">init-lambda</A><br>
+<A HREF="#run">init-step</A><br>
<A HREF="#run">integrator</A><br>
-<A HREF="#ld">ld_seed</A><br>
-<A HREF="#bond2">lincs_iter</A><br>
-<A HREF="#bond2">lincs_order</A><br>
-<A HREF="#bond2">lincs_warnangle</A><br>
+<A HREF="#ld">ld-seed</A><br>
+<A HREF="#bond2">lincs-iter</A><br>
+<A HREF="#bond2">lincs-order</A><br>
+<A HREF="#bond2">lincs-warnangle</A><br>
<A HREF="#bond2">morse</A><br>
<A HREF="#em">nbfgscorr</A><br>
<A HREF="#xmdrun">niter</A><br>
<A HREF="#out">nstvout</A><br>
<A HREF="#out">nstxout</A><br>
<A HREF="#out">nstxtcout</A><br>
-<A HREF="#nl">ns_type</A><br>
+<A HREF="#nl">ns-type</A><br>
<A HREF="#wall">nwall</A><br>
-<A HREF="#ewald">optimize_fft</A><br>
+<A HREF="#ewald">optimize-fft</A><br>
<A HREF="#nmr2">orire</A><br>
-<A HREF="#nmr2">orire_fc</A><br>
-<A HREF="#nmr2">orire_tau</A><br>
-<A HREF="#nmr2">orire_fitgrp</A><br>
+<A HREF="#nmr2">orire-fc</A><br>
+<A HREF="#nmr2">orire-tau</A><br>
+<A HREF="#nmr2">orire-fitgrp</A><br>
<A HREF="#nmr2">nstorireout</A><br>
<A HREF="#nl">pbc</A><br>
<A HREF="#pc">pcoupl</A><br>
<A HREF="#pc">pcoupltype</A><br>
-<A HREF="#nl">periodic_molecules</A><br>
-<A HREF="#ewald">pme_order</A><br>
+<A HREF="#nl">periodic-molecules</A><br>
+<A HREF="#ewald">pme-order</A><br>
<A HREF="#pull">pull</A><br>
-<A HREF="#pc">refcoord_scaling</A><br>
-<A HREF="#pc">ref_p</A><br>
-<A HREF="#tc">ref_t</A><br>
-<A HREF="#el2">rcoulomb_switch</A><br>
+<A HREF="#pc">refcoord-scaling</A><br>
+<A HREF="#pc">ref-p</A><br>
+<A HREF="#tc">ref-t</A><br>
+<A HREF="#el2">rcoulomb-switch</A><br>
<A HREF="#el2">rcoulomb</A><br>
<A HREF="#nl">rlist</A><br>
<A HREF="#nl">rlistlong</A><br>
<A HREF="#tpi">rtpi</A><br>
-<A HREF="#vdw">rvdw_switch</A><br>
+<A HREF="#vdw">rvdw-switch</A><br>
<A HREF="#vdw">rvdw</A><br>
-<A HREF="#free">sc_alpha</A><br>
-<A HREF="#free">sc_power</A><br>
-<A HREF="#free">sc_sigma</A><br>
-<A HREF="#bond2">shake_tol</A><br>
+<A HREF="#free">sc-alpha</A><br>
+<A HREF="#free">sc-power</A><br>
+<A HREF="#free">sc-sigma</A><br>
+<A HREF="#bond2">shake-tol</A><br>
<A HREF="#table">table-extension</A><br>
-<A HREF="#pc">tau_p</A><br>
-<A HREF="#tc">tau_t</A><br>
-<A HREF="#tc">tc_grps</A><br>
+<A HREF="#pc">tau-p</A><br>
+<A HREF="#tc">tau-t</A><br>
+<A HREF="#tc">tc-grps</A><br>
<A HREF="#tc">tcoupl</A><br>
<A HREF="#run">tinit</A><br>
<A HREF="#bond">continuation</A><br>
-<A HREF="#user">user1_grps</A><br>
-<A HREF="#user">user2_grps</A><br>
+<A HREF="#user">user1-grps</A><br>
+<A HREF="#user">user2-grps</A><br>
<A HREF="#user">userint1</A><br>
<A HREF="#user">userint2</A><br>
<A HREF="#user">userint3</A><br>
<A HREF="#user">userreal3</A><br>
<A HREF="#user">userreal4</A><br>
<A HREF="#el">vdwtype</A><br>
-<A HREF="#out">xtc_grps</A><br>
-<A HREF="#out">xtc_precision</A><br>
-<A HREF="#sa">zero_temp_time</A><br>
-<A HREF="#walls">wall_atomtype</A><br>
-<A HREF="#walls">wall_density</A><br>
-<A HREF="#walls">wall_ewald_zfac</A><br>
-<A HREF="#walls">wall_r_linpot</A><br>
-<A HREF="#walls">wall_type</A><br>
+<A HREF="#out">xtc-grps</A><br>
+<A HREF="#out">xtc-precision</A><br>
+<A HREF="#sa">zero-temp-time</A><br>
+<A HREF="#walls">wall-atomtype</A><br>
+<A HREF="#walls">wall-density</A><br>
+<A HREF="#walls">wall-ewald-zfac</A><br>
+<A HREF="#walls">wall-r-linpot</A><br>
+<A HREF="#walls">wall-type</A><br>
</multicol>
<hr>
#include "typedefs.h"
#include "smalloc.h"
#include "sysstuff.h"
-#include "errno.h"
+#include <errno.h>
#include "macros.h"
#include "string2.h"
#include "confio.h"
{
t_atoms *atoms;
gmx_bool bEnd;
- int nwanted,natoms,atnr,resnr,oldres,newres,shift;
+ int nwanted,natoms,atnr,resnr=0,oldres,newres,shift;
char anm[STRLEN],resnm[STRLEN];
char c1,c2;
double db1,db2,db3;
"Found more coordinates (%d) in %s than expected %d\n",
natoms,infile,nwanted);
if (atoms) {
- if (atoms && fr->bAtoms &&
+ if (fr->bAtoms &&
(sscanf(line,"%5d%c%5s%c%5s%7d",&resnr,&c1,resnm,&c2,anm,&atnr)
!= 6)) {
if (oldres>=0)
return natoms;
}
- int read_g96_conf(FILE *fp,const char *infile,t_trxframe *fr)
+ int read_g96_conf(FILE *fp,const char *infile,t_trxframe *fr, char *line)
{
t_symtab *symtab=NULL;
- char line[STRLEN+1];
gmx_bool bAtStart,bTime,bAtoms,bPos,bVel,bBox,bEnd,bFinished;
int natoms,nbp;
double db1,db2,db3,db4,db5,db6,db7,db8,db9;
gmx_fio_fclose (in);
}
- static gmx_bool get_w_conf(FILE *in,const char *infile,char *title,
- t_atoms *atoms, int *ndec, rvec x[],rvec *v, matrix box)
+ static gmx_bool get_w_conf(FILE *in, const char *infile, char *title,
+ t_symtab *symtab, t_atoms *atoms, int *ndec,
+ rvec x[], rvec *v, matrix box)
{
- t_symtab *symtab=NULL;
char name[6];
char line[STRLEN+1],*ptr;
char buf[256];
oldres = NOTSET; /* Unlikely number for the first residue! */
ddist = 0;
- if (!symtab) {
- snew(symtab,1);
- open_symtab(symtab);
- }
-
/* Read the title and number of atoms */
get_coordnum_fp(in,title,&natoms);
box[ZZ][XX] = y2;
box[ZZ][YY] = z2;
}
- close_symtab(symtab);
return bVel;
}
{
FILE *in;
int ndec;
-
+ t_symtab symtab;
+
/* open file */
in=gmx_fio_fopen(infile,"r");
- get_w_conf(in, infile, title, atoms, &ndec, x, v, box);
-
+ open_symtab(&symtab);
+ get_w_conf(in, infile, title, &symtab, atoms, &ndec, x, v, box);
+ /* We can't free the symbols, as they are still used in atoms, so
+ * the only choice is to leak them. */
+ free_symtab(&symtab);
+
gmx_fio_fclose(in);
}
gmx_bool gro_next_x_or_v(FILE *status,t_trxframe *fr)
{
t_atoms atoms;
+ t_symtab symtab;
char title[STRLEN],*p;
double tt;
- int ndec,i;
+ int ndec=0,i;
if (gmx_eof(status))
return FALSE;
+ open_symtab(&symtab);
atoms.nr=fr->natoms;
snew(atoms.atom,fr->natoms);
atoms.nres=fr->natoms;
snew(atoms.resinfo,fr->natoms);
snew(atoms.atomname,fr->natoms);
- fr->bV = get_w_conf(status,title,title,&atoms,&ndec,fr->x,fr->v,fr->box);
+ fr->bV = get_w_conf(status,title,title,&symtab,&atoms,&ndec,fr->x,fr->v,fr->box);
fr->bPrec = TRUE;
fr->prec = 1;
/* prec = 10^ndec: */
sfree(atoms.atom);
sfree(atoms.resinfo);
sfree(atoms.atomname);
+ done_symtab(&symtab);
if ((p=strstr(title,"t=")) != NULL) {
p+=2;
FILE *in;
int ftp,tpxver,tpxgen;
t_trxframe fr;
+ char g96_line[STRLEN+1];
ftp=fn2ftp(infile);
range_check(ftp,0,efNR);
fr.x = NULL;
fr.v = NULL;
fr.f = NULL;
- *natoms=read_g96_conf(in,infile,&fr);
+ *natoms=read_g96_conf(in,infile,&fr,g96_line);
gmx_fio_fclose(in);
break;
case efXYZ:
t_trxframe fr;
int i,ftp,natoms;
real d;
+ char g96_line[STRLEN+1];
if (atoms->nr == 0)
fprintf(stderr,"Warning: Number of atoms in %s is 0\n",infile);
fr.v = v;
fr.f = NULL;
in = gmx_fio_fopen(infile,"r");
- read_g96_conf(in, infile, &fr);
+ read_g96_conf(in, infile, &fr, g96_line);
gmx_fio_fclose(in);
copy_mat(fr.box,box);
break;
* structure, so we should not free them. But there is no place to put the
* symbols; the only choice is to leak the memory...
* So we clear the symbol table before freeing the topology structure. */
- open_symtab(&top.symtab);
+ free_symtab(&top.symtab);
done_top(&top);
break;
#include <config.h>
#endif
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
#include <thread_mpi.h>
#endif
#include "smalloc.h"
#include "string2.h"
#include "macros.h"
-#include "time.h"
+#include <time.h>
#include "random.h"
#include "statutil.h"
#include "copyrite.h"
"Phys. Chem. Chem. Phys.",
13, 2011, "169-181" },
{ "Caleman2011b",
- "C. Caleman and M. Hong and J. S. Hub and L. T. da Costa and P. J. van Maaren and D. van der Spoel",
- "Force Field Benchmark 1: Density, Heat of Vaporization, Heat Capacity, Surface Tension and Dielectric Constant of 152 Organic Liquids",
- "Submitted",
- 0, 2011, "" },
+ "C. Caleman and P. J. van Maaren and M. Hong and J. S. Hub and L. T. da Costa and D. van der Spoel",
+ "Force Field Benchmark of Organic Liquids: Density, Enthalpy of Vaporization, Heat Capacities, Surface Tension, Isothermal Compressibility, Volumetric Expansion Coefficient, and Dielectric Constant",
+ "J. Chem. Theo. Comp.",
+ 8, 2012, "61" },
{ "Lindahl2001a",
"E. Lindahl and B. Hess and D. van der Spoel",
"GROMACS 3.0: A package for molecular simulation and trajectory analysis",
"O. Engin, A. Villa, M. Sayar and B. Hess",
"Driving Forces for Adsorption of Amphiphilic Peptides to Air-Water Interface",
"J. Phys. Chem. B",
- 0, 2010, "???" },
+ 114, 2010, "11093" },
+ { "Fritsch12",
+ "S. Fritsch, C. Junghans and K. Kremer",
+ "Adaptive molecular simulation study on structure formation of toluene around C60 using Gromacs",
+ "J. Chem. Theo. Comp.",
+ 8, 2012, "398" },
+ { "Junghans10",
+ "C. Junghans and S. Poblete",
+ "A reference implementation of the adaptive resolution scheme in ESPResSo",
+ "Comp. Phys. Comm.",
+ 181, 2010, "1449" },
{ "Wang2010",
"H. Wang, F. Dommert, C.Holm",
"Optimizing working parameters of the smooth particle mesh Ewald algorithm in terms of accuracy and efficiency",
"J. Chem. Phys. B",
- 133, 2010, "034117"
- },
+ 133, 2010, "034117" },
+ { "Sugita1999a",
+ "Y. Sugita, Y. Okamoto",
+ "Replica-exchange molecular dynamics method for protein folding",
+ "Chem. Phys. Lett.",
+ 314, 1999, "141-151" },
+ { "Kutzner2011",
+ "C. Kutzner and J. Czub and H. Grubmuller",
+ "Keep it Flexible: Driving Macromolecular Rotary Motions in Atomistic Simulations with GROMACS",
+ "J. Chem. Theory Comput.",
+ 7, 2011, "1381-1393" },
+ { "Hoefling2011",
+ "M. Hoefling, N. Lima, D. Haenni, C.A.M. Seidel, B. Schuler, H. Grubmuller",
+ "Structural Heterogeneity and Quantitative FRET Efficiency Distributions of Polyprolines through a Hybrid Atomistic Simulation and Monte Carlo Approach",
+ "PLoS ONE",
+ 6, 2011, "e19791" },
+ { "Hockney1988",
+ "R. W. Hockney and J. W. Eastwood",
+ "Computer simulation using particles",
+ "IOP, Bristol",
+ 1, 1988, "1" },
+ { "Ballenegger2012",
+ "V. Ballenegger, J.J. Cerda, and C. Holm",
+ "How to Convert SPME to P3M: Influence Functions and Error Estimates",
+ "J. Chem. Theory Comput.",
+ 8, 2012, "936-947" },
+ { "Garmay2012",
+ "Garmay Yu, Shvetsov A, Karelov D, Lebedev D, Radulescu A, Petukhov M, Isaev-Ivanov V",
+ "Correlated motion of protein subdomains and large-scale conformational flexibility of RecA protein filament",
+ "Journal of Physics: Conference Series",
+ 340, 2012, "012094" }
};
#define NSTR (int)asize(citedb)
fprintf(fp, "Precision: single\n");
#endif
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
fprintf(fp, "Parallellization: thread_mpi\n");
#elif defined(GMX_MPI)
fprintf(fp, "Parallellization: MPI\n");
#ifdef GMX_FFT_FFTPACK
fprintf(fp, "FFT Library: fftpack\n");
-#elif defined(GMX_FFT_FFTW2)
- fprintf(fp, "FFT Library: fftw2\n");
#elif defined(GMX_FFT_FFTW3)
fprintf(fp, "FFT Library: fftw3\n");
#elif defined(GMX_FFT_MKL)
#include <ctype.h>
#include <string.h>
+#include <assert.h>
#include "sysstuff.h"
#include "strdb.h"
#include "futil.h"
for(i=0;i<rt->n;i++)
{
- free(rt->resname[i]);
- free(rt->restype[i]);
+ sfree(rt->resname[i]);
+ sfree(rt->restype[i]);
}
- free(rt);
+ sfree(rt->resname);
+ sfree(rt->restype);
+ sfree(rt);
return 0;
}
void analyse(t_atoms *atoms,t_blocka *gb,char ***gn,gmx_bool bASK,gmx_bool bVerb)
{
- gmx_residuetype_t rt;
+ gmx_residuetype_t rt=NULL;
char *resnm;
atom_id *aid;
const char ** restype;
/* For every residue, get a pointer to the residue type name */
gmx_residuetype_init(&rt);
+ assert(rt);
snew(restype,atoms->nres);
ntypes = 0;
b->index[b->nr]=b->index[b->nr-1];
(*grpname)[b->nr-1]=strdup(str);
} else {
+ if (b->nr==0)
+ {
+ gmx_fatal(FARGS,"The first header of your indexfile is invalid");
+ }
pt=line;
- while ((i=sscanf(pt,"%s",str)) == 1) {
+ while (sscanf(pt,"%s",str) == 1) {
i=b->index[b->nr];
if (i>=maxentries) {
maxentries+=1024;
fprintf(out,"CONECT%5d%5d\n",gc->conect[i].ai+1,gc->conect[i].aj+1);
}
}
+
+ gmx_residuetype_destroy(rt);
}
void write_pdbfile(FILE *out,const char *title, t_atoms *atoms,rvec x[],
for(k=0; (k<4); k++,j++) anm[k]=line[j];
anm[k]=nc;
strcpy(anm_copy,anm);
+ rtrim(anm_copy);
atomnumber = NOTSET;
trim(anm);
altloc=line[j];
{
c=line+6;
/* skip HEADER or TITLE and spaces */
- while (c && (c[0]!=' ')) c++;
- while (c && (c[0]==' ')) c++;
+ while (c[0]!=' ') c++;
+ while (c[0]==' ') c++;
/* truncate after title */
d=strstr(c," ");
if (d)
c=line;
}
/* skip 'MOLECULE:' and spaces */
- while (c && (c[0]!=' ')) c++;
- while (c && (c[0]==' ')) c++;
+ while (c[0]!=' ') c++;
+ while (c[0]==' ') c++;
/* truncate after title */
d=strstr(c," ");
if (d)
#include <smalloc.h>
#include <string2.h>
#include <vec.h>
+#include <assert.h>
#include <indexutil.h>
#include <poscalc.h>
*
* If called more than once, memory is (re)allocated to ensure that the
* maximum of the \p isize values can be stored.
+ *
+ * Allocation of POS_VALUE selection elements is a special case, and is
+ * handled by alloc_selection_pos_data().
*/
static gmx_bool
alloc_selection_data(t_selelem *sel, int isize, gmx_bool bChildEval)
{
child = child->child;
}
- nalloc = (sel->v.type == POS_VALUE) ? child->v.u.p->nr : child->v.nr;
- }
- /* For positions, we actually want to allocate just a single structure
- * for nalloc positions. */
- if (sel->v.type == POS_VALUE)
- {
- isize = nalloc;
- nalloc = 1;
+ nalloc = child->v.nr;
}
/* Allocate memory for sel->v.u if needed */
if (sel->flags & SEL_ALLOCVAL)
{
_gmx_selvalue_reserve(&sel->v, nalloc);
}
- /* Reserve memory inside group and position structures if
- * SEL_ALLOCDATA is set. */
- if (sel->flags & SEL_ALLOCDATA)
+ /* Reserve memory inside group structure if SEL_ALLOCDATA is set. */
+ if ((sel->flags & SEL_ALLOCDATA) && sel->v.type == GROUP_VALUE)
{
- if (sel->v.type == GROUP_VALUE)
- {
- gmx_ana_index_reserve(sel->v.u.g, isize);
- }
- else if (sel->v.type == POS_VALUE)
- {
- gmx_ana_pos_reserve(sel->v.u.p, isize, 0);
- }
+ gmx_ana_index_reserve(sel->v.u.g, isize);
}
return TRUE;
}
+ /*! \brief
+ * Allocates memory for storing the evaluated value of a selection element.
+ *
+ * \param sel Selection element to initialize.
+ *
+ * Allocation of POS_VALUE selection elements is a special case, and is
+ * handled by this function instead of by alloc_selection_data().
+ */
+ static void
+ alloc_selection_pos_data(t_selelem *sel)
+ {
+ t_selelem *child;
+ int nalloc, isize;
+
+ if (sel->mempool)
+ {
+ return;
+ }
+
+ child = sel;
+ if (sel->type == SEL_SUBEXPRREF)
+ {
+ child = sel->child->child;
+ }
+ nalloc = child->v.u.p->nr;
+ isize = child->v.u.p->m.b.nra;
+
+ /* For positions, we want to allocate just a single structure
+ * for nalloc positions. */
+ if (sel->flags & SEL_ALLOCVAL)
+ {
+ _gmx_selvalue_reserve(&sel->v, 1);
+ }
+ sel->v.nr = 1;
+ /* Reserve memory inside position structure if SEL_ALLOCDATA is set. */
+ if (sel->flags & SEL_ALLOCDATA)
+ {
+ gmx_ana_pos_reserve(sel->v.u.p, nalloc, isize);
+ }
+ }
+
/*! \brief
* Replace the evaluation function of each element in the subtree.
*
static void
init_item_evaloutput(t_selelem *sel)
{
+ assert(!(sel->child == NULL &&
+ (sel->type == SEL_SUBEXPRREF || sel->type == SEL_SUBEXPR)));
+
/* Process children. */
if (sel->type != SEL_SUBEXPRREF)
{
&& ((sel->cdata->flags & SEL_CDATA_SIMPLESUBEXPR)
|| (sel->cdata->flags & SEL_CDATA_FULLEVAL)))
{
+ assert(sel->child);
sel->cdata->gmin = sel->child->cdata->gmin;
sel->cdata->gmax = sel->child->cdata->gmax;
}
* \param[in,out] sc Selection collection data.
*
* The evaluation group of each \ref SEL_ROOT element corresponding to a
- * selection in \p sc is set to \p gall. The same is done for \ref SEL_ROOT
- * elements corresponding to subexpressions that need full evaluation.
+ * selection in \p sc is set to NULL. The evaluation grop for \ref SEL_ROOT
+ * elements corresponding to subexpressions that need full evaluation is set
+ * to \c sc->gall.
*/
static void
initialize_evalgrps(gmx_ana_selcollection_t *sc)
while (root)
{
if (root->child->type != SEL_SUBEXPR
- || (root->child->cdata->flags & SEL_CDATA_FULLEVAL))
+ || (root->child->v.type != GROUP_VALUE && !(root->flags & SEL_ATOMVAL)))
+ {
+ gmx_ana_index_set(&root->u.cgrp, -1, 0, root->u.cgrp.name, 0);
+ }
+ else if (root->child->cdata->flags & SEL_CDATA_FULLEVAL)
{
gmx_ana_index_set(&root->u.cgrp, sc->gall.isize, sc->gall.index,
root->u.cgrp.name, 0);
{
return rc;
}
- if (sel->v.type != POS_VALUE && sel->v.type != GROUP_VALUE)
+ if (sel->v.type != POS_VALUE && sel->v.type != GROUP_VALUE
+ && !(sel->flags & SEL_VARNUMVAL))
{
alloc_selection_data(sel, isize, TRUE);
}
case SEL_EXPRESSION:
case SEL_MODIFIER:
+ assert(g);
rc = _gmx_sel_evaluate_method_params(data, sel, g);
if (rc != 0)
{
return rc;
}
- rc = init_method(sel, data->top, g->isize);
+ rc = init_method(sel, data->top, g ? g->isize : 0);
if (rc != 0)
{
return rc;
{
rc = sel->cdata->evaluate(data, sel, g);
}
- if (bDoMinMax)
+ if (bDoMinMax && g)
{
gmx_ana_index_copy(sel->cdata->gmax, g, TRUE);
}
}
else if (sel->u.cgrp.isize == 0)
{
+ assert(g);
gmx_ana_index_reserve(&sel->u.cgrp, g->isize);
rc = sel->cdata->evaluate(data, sel, g);
if (bDoMinMax)
/* The subexpression should have been evaluated if g is NULL
* (i.e., this is a method parameter or a direct value of a
* selection). */
- alloc_selection_data(sel, sel->child->cdata->gmax->isize, TRUE);
+ if (sel->v.type == POS_VALUE)
+ {
+ alloc_selection_pos_data(sel);
+ }
+ else
+ {
+ alloc_selection_data(sel, sel->child->cdata->gmax->isize, TRUE);
+ }
}
rc = sel->cdata->evaluate(data, sel, g);
if (rc != 0)
static void
postprocess_item_subexpressions(t_selelem *sel)
{
+ assert(!(sel->child == NULL &&
+ (sel->type == SEL_SUBEXPRREF || sel->type == SEL_SUBEXPR)));
+
/* Process children. */
if (sel->type != SEL_SUBEXPRREF)
{
sel->u.cgrp.name = name;
sel->evaluate = &_gmx_sel_evaluate_subexpr_staticeval;
- if (sel->cdata)
- {
- sel->cdata->evaluate = sel->evaluate;
- }
+ sel->cdata->evaluate = sel->evaluate;
+
_gmx_selelem_free_values(sel->child);
sel->child->mempool = NULL;
_gmx_selvalue_setstore(&sel->child->v, sel->v.u.ptr);
#include <maths.h>
#include <smalloc.h>
#include <vec.h>
+#include <assert.h>
#include <indexutil.h>
#include <poscalc.h>
return rc;
}
sel->v.nr = sel->child->v.nr;
- gmx_ana_index_set(&sel->u.cgrp, g->isize, g->index, sel->u.cgrp.name, 0);
+ if (!g)
+ {
+ sel->u.cgrp.isize = -1;
+ }
+ else
+ {
+ gmx_ana_index_set(&sel->u.cgrp, g->isize, g->index, sel->u.cgrp.name, 0);
+ }
}
return 0;
}
{
_gmx_sel_mempool_free_group(data->mp, &gmiss);
}
+ gmiss.name = NULL;
}
if (gmiss.isize > 0)
{
{
if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
{
- sel->v.u.i[k] = sel->v.u.i[j--];
+ sel->v.u.i[k] = sel->child->v.u.i[j--];
}
else
{
- sel->v.u.i[k] = sel->child->v.u.i[i--];
+ sel->v.u.i[k] = sel->v.u.i[i--];
}
}
break;
{
if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
{
- sel->v.u.r[k] = sel->v.u.r[j--];
+ sel->v.u.r[k] = sel->child->v.u.r[j--];
}
else
{
- sel->v.u.r[k] = sel->child->v.u.r[i--];
+ sel->v.u.r[k] = sel->v.u.r[i--];
}
}
break;
{
if (i < 0 || (j >= 0 && sel->u.cgrp.index[i] < gmiss.index[j]))
{
- sel->v.u.s[k] = sel->v.u.s[j--];
+ sel->v.u.s[k] = sel->child->v.u.s[j--];
}
else
{
- sel->v.u.s[k] = sel->child->v.u.s[i--];
+ sel->v.u.s[k] = sel->v.u.s[i--];
}
}
break;
_gmx_sel_evaluate_arithmetic(gmx_sel_evaluate_t *data, t_selelem *sel,
gmx_ana_index_t *g)
{
- t_selelem *left, *right;
int n, i, i1, i2;
real lval, rval=0., val=0.;
int rc;
+ gmx_bool bArithNeg;
- left = sel->child;
- right = left->next;
+ t_selelem *const left = sel->child;
+ t_selelem *const right = left->next;
if (left->mempool)
{
n = (sel->flags & SEL_SINGLEVAL) ? 1 : g->isize;
sel->v.nr = n;
+
+ bArithNeg = (sel->u.arith.type == ARITH_NEG);
+ assert(right || bArithNeg);
for (i = i1 = i2 = 0; i < n; ++i)
{
lval = left->v.u.r[i1];
- if (sel->u.arith.type != ARITH_NEG)
+ if (!bArithNeg)
{
rval = right->v.u.r[i2];
}
{
++i1;
}
- if (sel->u.arith.type != ARITH_NEG && !(right->flags & SEL_SINGLEVAL))
+ if (!bArithNeg && !(right->flags & SEL_SINGLEVAL))
{
++i2;
}
extern gmx_ana_selmethod_t sm_resindex;
extern gmx_ana_selmethod_t sm_molindex;
extern gmx_ana_selmethod_t sm_atomname;
+ extern gmx_ana_selmethod_t sm_pdbatomname;
extern gmx_ana_selmethod_t sm_atomtype;
extern gmx_ana_selmethod_t sm_resname;
extern gmx_ana_selmethod_t sm_insertcode;
{"mol", &sm_molindex},
{"molecule", &sm_molindex},
{NULL, &sm_atomname},
+ {"name", &sm_atomname},
+ {NULL, &sm_pdbatomname},
+ {"pdbname", &sm_pdbatomname},
{NULL, &sm_atomtype},
+ {"type", &sm_atomtype},
{NULL, &sm_resname},
{NULL, &sm_insertcode},
{NULL, &sm_chain},
if (nparams > 0 && !param)
{
report_error(fp, name, "error: missing parameter data");
- bOk = FALSE;
return FALSE;
}
if (nparams == 0 && param)
report_error(fp, method->name, "error: outinit should be provided because the method has POS_VALUE");
bOk = FALSE;
}
+ /* Check presence of outinit for variable output count methods */
+ if ((method->flags & SMETH_VARNUMVAL) && !method->outinit)
+ {
+ report_error(fp, method->name, "error: outinit should be provided because the method has SMETH_VARNUMVAL");
+ bOk = FALSE;
+ }
/* Warn of dynamic callbacks in static methods */
if (!(method->flags & SMETH_MODIFIER))
{
* \f[\sin^2 \frac{\Delta \phi}{2} = \frac{\sin^2 \frac{\alpha}{2} - \sin^2 \frac{\theta - \theta'}{2}}{\sin \theta \sin \theta'}\f]
* (distance in spherical geometry),
* where \f$\theta'\f$ is the zenith angle of the bin edge.
+ * Treat zenith angle bins that are completely covered by the cone (in the
+ * case that the cone is centered close to the pole) as a special case.
* -# Using the values calculated above, loop through the azimuthal bins that
* are partially or completely covered by the cone and update them.
*
/** Initializes the \p insolidangle selection method. */
static int
init_insolidangle(t_topology *top, int npar, gmx_ana_selparam_t *param, void *data);
-/** Sets the COM/COG data for the \p insolidangle selection method. */
-static void
-set_comg_insolidangle(gmx_ana_pos_t *pos, void *data);
/** Frees the data allocated for the \p insolidangle selection method. */
static void
free_data_insolidangle(void *data);
rvec x)
{
real pdelta, phi1, phi2;
- int pbin1, pbin2, pbin;
+ int pbin1, pbin2, pbiniter, pbin;
/* Find the edges of the bins affected */
pdelta = max(max(pdelta1, pdelta2), pdeltamax);
phi1 = phi - pdelta;
- if (phi1 < -M_PI)
+ if (phi1 >= -M_PI)
{
- phi1 += M_2PI;
+ pbin = find_partition_bin(&surf->tbin[tbin], phi1);
+ pbin1 = pbin;
+ }
+ else
+ {
+ pbin = find_partition_bin(&surf->tbin[tbin], phi1 + M_2PI);
+ pbin1 = pbin - surf->tbin[tbin].n;
}
phi2 = phi + pdelta;
- if (phi2 > M_PI)
+ if (phi2 <= M_PI)
+ {
+ pbin2 = find_partition_bin(&surf->tbin[tbin], phi2);
+ }
+ else
+ {
+ pbin2 = find_partition_bin(&surf->tbin[tbin], phi2 - M_2PI);
+ pbin2 += surf->tbin[tbin].n;
+ }
+ ++pbin2;
+ if (pbin2 - pbin1 > surf->tbin[tbin].n)
{
- phi2 -= M_2PI;
+ pbin2 = pbin1 + surf->tbin[tbin].n;
}
- pbin1 = find_partition_bin(&surf->tbin[tbin], phi1);
- pbin2 = find_partition_bin(&surf->tbin[tbin], phi2);
/* Find the edges of completely covered region */
pdelta = min(pdelta1, pdelta2);
phi1 = phi - pdelta;
}
phi2 = phi + pdelta;
/* Loop over all affected bins */
- pbin = pbin1;
- do
+ for (pbiniter = pbin1; pbiniter != pbin2; ++pbiniter, ++pbin)
{
/* Wrap bin around if end reached */
if (pbin == surf->tbin[tbin].n)
add_surface_point(surf, tbin, pbin, x);
}
}
- while (pbin++ != pbin2); /* Loop including pbin2 */
}
/*!
theta2 = (tbin+1) * surf->tbinsize;
if (theta2 > theta + surf->angcut)
{
+ /* The circle is completely outside the cone */
pdelta2 = 0;
}
- else if (tbin == surf->ntbins - 1)
+ else if (theta2 <= -(theta - surf->angcut)
+ || theta2 >= M_2PI - (theta + surf->angcut)
+ || tbin == surf->ntbins - 1)
{
+ /* The circle is completely inside the cone, or we are in the
+ * 360 degree bin covering the pole. */
pdelta2 = M_PI;
}
else
{
+ /* TODO: This formula is numerically unstable if theta is very
+ * close to the pole. In practice, it probably does not matter
+ * much, but it would be nicer to adjust the theta bin boundaries
+ * such that the case above catches this instead of falling through
+ * here. */
pdelta2 = 2*asin(sqrt(
(sqr(sin(surf->angcut/2)) - sqr(sin((theta2-theta)/2))) /
(sin(theta) * sin(theta2))));
#include <ctype.h>
+#include <assert.h>
#include "copyrite.h"
#include "sysstuff.h"
#include "macros.h"
#include "mtop_util.h"
#include "gmxfio.h"
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
#include "thread_mpi.h"
#endif
static const char *program_name=NULL;
static char *cmd_line=NULL;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
/* For now, some things here are simply not re-entrant, so
we have to actively lock them. */
static tMPI_Thread_mutex_t init_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
const char *ShortProgram(void)
{
const char *pr,*ret;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_lock(&init_mutex);
#endif
pr=ret=program_name;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_unlock(&init_mutex);
#endif
if ((pr=strrchr(ret,DIR_SEPARATOR)) != NULL)
ret=pr+1;
- /* Strip away the libtool prefix if it's still there. */
- if(strlen(ret) > 3 && !strncmp(ret, "lt-", 3))
- ret = ret + 3;
return ret;
}
const char *Program(void)
{
const char *ret;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_lock(&init_mutex);
#endif
ret=program_name;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_unlock(&init_mutex);
#endif
return ret;
const char *command_line(void)
{
const char *ret;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_lock(&init_mutex);
#endif
ret=cmd_line;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_unlock(&init_mutex);
#endif
return ret;
void set_program_name(const char *argvzero)
{
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_lock(&init_mutex);
#endif
- /* When you run a dynamically linked program before installing
- * it, libtool uses wrapper scripts and prefixes the name with "lt-".
- * Until libtool is fixed to set argv[0] right, rip away the prefix:
- */
if (program_name == NULL)
{
- if(strlen(argvzero)>3 && !strncmp(argvzero,"lt-",3))
- program_name=strdup(argvzero+3);
- else
- program_name=strdup(argvzero);
+ /* if filename has file ending (e.g. .exe) then strip away */
+ char* extpos=strrchr(argvzero,'.');
+ if(extpos > strrchr(argvzero,DIR_SEPARATOR))
+ {
+ program_name=gmx_strndup(argvzero,extpos-argvzero);
+ }
+ else
+ {
+ program_name=gmx_strdup(argvzero);
+ }
}
if (program_name == NULL)
program_name="GROMACS";
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_unlock(&init_mutex);
#endif
}
int i;
size_t cmdlength;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_lock(&init_mutex);
#endif
if (cmd_line==NULL)
strcat(cmd_line," ");
}
}
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_unlock(&init_mutex);
#endif
static void set_default_time_unit(const char *time_list[], gmx_bool bCanTime)
{
- int i,j;
+ int i=0,j;
- const char *select;
+ const char *select = NULL;
if (bCanTime)
{
}
}
}
- if (!bCanTime || select == NULL || strcmp(time_list[i], select) != 0)
+ if (!bCanTime || select == NULL ||
+ time_list[i]==NULL || strcmp(time_list[i], select) != 0)
{
/* Set it to the default: ps */
i = 1;
static void usage(const char *type,const char *arg)
{
- if (arg != NULL)
- gmx_fatal(FARGS,"Expected %s argument for option %s\n",type,arg);
+ assert(arg);
+ gmx_fatal(FARGS,"Expected %s argument for option %s\n",type,arg);
}
int iscan(int argc,char *argv[],int *i)
/* The some system, e.g. the catamount kernel on cray xt3 do not have nice(2). */
if (nicelevel != 0 && !bExit)
{
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
static gmx_bool nice_set=FALSE; /* only set it once */
tMPI_Thread_mutex_lock(&init_mutex);
if (!nice_set)
{
#endif
i=nice(nicelevel); /* assign ret value to avoid warnings */
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
nice_set=TRUE;
}
tMPI_Thread_mutex_unlock(&init_mutex);
gmx_cmd(argv[1]);
}
}
- if (bExit) {
- if (gmx_parallel_env_initialized())
- /*gmx_abort(gmx_node_rank(),gmx_node_num(),0);*/
- gmx_finalize();
+ if (bExit)
+ {
+ gmx_finalize_par();
+
exit(0);
}
#undef FF
#include <ctype.h>
#include "sysstuff.h"
+#include "typedefs.h"
+#include "vmdio.h"
#include "string2.h"
#include "smalloc.h"
#include "pbc.h"
#include "confio.h"
#include "checkpoint.h"
#include "wgms.h"
-#include "vmdio.h"
#include <math.h>
/* defines for frame counter output */
int NATOMS;
double DT,BOX[3];
gmx_bool bReadBox;
+ char *persistent_line; /* Persistent line for reading g96 trajectories */
};
static void initcount(t_trxstatus *status)
status->xframe=NULL;
status->fio=NULL;
status->__frame=-1;
+ status->persistent_line=NULL;
}
for(i=0; i<nind; i++)
copy_rvec(fr->f[ind[i]],fout[i]);
}
+ /* no break */
case efXTC:
case efG87:
if (fr->bX) {
case efTRR:
if (vout) sfree(vout);
if (fout) sfree(fout);
+ /* no break */
case efXTC:
case efG87:
sfree(xout);
fr->bTime=TRUE;
fr->time=sh.t;
fr->bLambda = TRUE;
+ fr->bFepState = TRUE;
fr->lambda = sh.lambda;
fr->bBox = sh.box_size>0;
if (fr->flags & (TRX_READ_X | TRX_NEED_X)) {
return bRet;
}
-static void choose_ff(FILE *fp)
+static void choose_file_format(FILE *fp)
{
int i,m,c;
int rc;
initcount(status);
clear_mat(box);
- choose_ff(fp);
+ choose_file_format(fp);
for(m=0; (m<DIM); m++)
box[m][m]=status->BOX[m];
/* Checkpoint files can not contain mulitple frames */
break;
case efG96:
- gmx_fatal(FARGS,
- "Reading trajectories in .g96 format is broken. Please use\n"
- "a different file format.");
- read_g96_conf(gmx_fio_getfp(status->fio),NULL,fr);
+ read_g96_conf(gmx_fio_getfp(status->fio),NULL,fr,
+ status->persistent_line);
bRet = (fr->natoms > 0);
break;
case efG87:
bRet = gro_next_x_or_v(gmx_fio_getfp(status->fio),fr);
break;
default:
-#ifdef GMX_DLOPEN
+#ifdef GMX_USE_PLUGINS
bRet = read_next_vmd_frame(dummy,fr);
#else
gmx_fatal(FARGS,"DEATH HORROR in read_next_frame ftp=%s,status=%s",
}
if (bRet) {
- bMissingData = ((fr->flags & TRX_NEED_X && !fr->bX) ||
- (fr->flags & TRX_NEED_V && !fr->bV) ||
- (fr->flags & TRX_NEED_F && !fr->bF));
+ bMissingData = (((fr->flags & TRX_NEED_X) && !fr->bX) ||
+ ((fr->flags & TRX_NEED_V) && !fr->bV) ||
+ ((fr->flags & TRX_NEED_F) && !fr->bF));
bSkip = FALSE;
if (!bMissingData) {
ct=check_times2(fr->time,fr->t0,fr->tpf,fr->tppf,fr->bDouble);
- if (ct == 0 || (fr->flags & TRX_DONT_SKIP && ct<0)) {
+ if (ct == 0 || ((fr->flags & TRX_DONT_SKIP) && ct<0)) {
printcount(status, oenv,fr->time,FALSE);
} else if (ct > 0)
bRet = FALSE;
break;
case efG96:
/* Can not rewind a compressed file, so open it twice */
- read_g96_conf(gmx_fio_getfp(fio),fn,fr);
+ if (!(*status)->persistent_line)
+ {
+ /* allocate the persistent line */
+ snew((*status)->persistent_line, STRLEN+1);
+ }
+ read_g96_conf(gmx_fio_getfp(fio),fn,fr, (*status)->persistent_line);
gmx_fio_close(fio);
clear_trxframe(fr,FALSE);
if (flags & (TRX_READ_X | TRX_NEED_X))
bFirst = FALSE;
break;
default:
-#ifdef GMX_DLOPEN
+#ifdef GMX_USE_PLUGINS
fprintf(stderr,"The file format of %s is not a known trajectory format to GROMACS.\n"
"Please make sure that the file is a trajectory!\n"
"GROMACS will now assume it to be a trajectory and will try to open it using the VMD plug-ins.\n"
}
#else
gmx_fatal(FARGS,"Not supported in read_first_frame: %s. Please make sure that the file is a trajectory.\n"
- "GROMACS is not compiled with DLOPEN. Thus it cannot read non-GROMACS trajectory formats.\n"
- "Please compile with DLOPEN support if you want to read non-GROMACS trajectory formats.\n",fn);
+ "GROMACS is not compiled with plug-in support. Thus it cannot read non-GROMACS trajectory formats using the VMD plug-ins.\n"
+ "Please compile with plug-in support if you want to read non-GROMACS trajectory formats.\n",fn);
#endif
break;
}
clear_rvec(fr->v[i]);
}
-int read_first_v(const output_env_t oenv, t_trxstatus **status,const char *fn,
- real *t, rvec **v,matrix box)
-{
- t_trxframe fr;
-
- read_first_frame(oenv,status,fn,&fr,TRX_NEED_V);
- *t = fr.time;
- clear_v(&fr);
- *v = fr.v;
- copy_mat(fr.box,box);
-
- return fr.natoms;
-}
-
-gmx_bool read_next_v(const output_env_t oenv,t_trxstatus *status,real *t,
- int natoms,rvec v[], matrix box)
-{
- t_trxframe fr;
- gmx_bool bRet;
-
- clear_trxframe(&fr,TRUE);
- fr.flags = TRX_NEED_V;
- fr.natoms = natoms;
- fr.time = *t;
- fr.v = v;
- bRet = read_next_frame(oenv,status,&fr);
- *t = fr.time;
- clear_v(&fr);
- copy_mat(fr.box,box);
-
- return bRet;
-}
#include "pbc.h"
#include <string.h>
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
#include "thread_mpi.h"
#endif
static gmx_bool bOverAllocDD=FALSE;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
static tMPI_Thread_mutex_t over_alloc_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
#endif
void set_over_alloc_dd(gmx_bool set)
{
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_lock(&over_alloc_mutex);
/* we just make sure that we don't set this at the same time.
We don't worry too much about reading this rarely-set variable */
#endif
bOverAllocDD = set;
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
tMPI_Thread_mutex_unlock(&over_alloc_mutex);
#endif
}
void init_inputrec(t_inputrec *ir)
{
- memset(ir,0,(size_t)sizeof(*ir));
+ memset(ir,0,(size_t)sizeof(*ir));
+ snew(ir->fepvals,1);
+ snew(ir->expandedvals,1);
+ snew(ir->simtempvals,1);
}
void stupid_fill_block(t_block *grp,int natom,gmx_bool bOneIndexGroup)
sfree(at->atomname);
sfree(at->atomtype);
sfree(at->atomtypeB);
+ if (at->pdbinfo)
+ sfree(at->pdbinfo);
}
void done_atomtypes(t_atomtypes *atype)
}
}
-static void init_ekinstate(ekinstate_t *eks)
+static void zero_ekinstate(ekinstate_t *eks)
{
eks->ekin_n = 0;
eks->ekinh = NULL;
for(i=0; i<state->ngtc; i++)
{
for (j=0;j<state->nhchainlength;j++)
- {
+ {
state->nosehoover_xi[i*state->nhchainlength + j] = 0.0;
state->nosehoover_vxi[i*state->nhchainlength + j] = 0.0;
}
}
-void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength)
+void init_state(t_state *state, int natoms, int ngtc, int nnhpres, int nhchainlength, int nlambda)
{
int i;
state->nrng = 0;
state->flags = 0;
state->lambda = 0;
+ snew(state->lambda,efptNR);
+ for (i=0;i<efptNR;i++)
+ {
+ state->lambda[i] = 0;
+ }
state->veta = 0;
clear_mat(state->box);
clear_mat(state->box_rel);
state->sd_X = NULL;
state->cg_p = NULL;
- init_ekinstate(&state->ekinstate);
+ zero_ekinstate(&state->ekinstate);
init_energyhistory(&state->enerhist);
+ init_df_history(&state->dfhist,nlambda,0);
+
state->ddp_count = 0;
state->ddp_count_cg_gl = 0;
state->cg_gl = NULL;
return max(cutoff1,cutoff2);
}
}
+
+extern void init_df_history(df_history_t *dfhist, int nlambda, real wl_delta)
+{
+ int i;
+
+ dfhist->bEquil = 0;
+ dfhist->nlambda = nlambda;
+ dfhist->wl_delta = wl_delta;
+ snew(dfhist->sum_weights,dfhist->nlambda);
+ snew(dfhist->sum_dg,dfhist->nlambda);
+ snew(dfhist->sum_minvar,dfhist->nlambda);
+ snew(dfhist->sum_variance,dfhist->nlambda);
+ snew(dfhist->n_at_lam,dfhist->nlambda);
+ snew(dfhist->wl_histo,dfhist->nlambda);
+
+ /* allocate transition matrices here */
+ snew(dfhist->Tij,dfhist->nlambda);
+ snew(dfhist->Tij_empirical,dfhist->nlambda);
+
+ for (i=0;i<dfhist->nlambda;i++) {
+ snew(dfhist->Tij[i],dfhist->nlambda);
+ snew(dfhist->Tij_empirical[i],dfhist->nlambda);
+ }
+
+ snew(dfhist->accum_p,dfhist->nlambda);
+ snew(dfhist->accum_m,dfhist->nlambda);
+ snew(dfhist->accum_p2,dfhist->nlambda);
+ snew(dfhist->accum_m2,dfhist->nlambda);
+
+ for (i=0;i<dfhist->nlambda;i++) {
+ snew((dfhist->accum_p)[i],dfhist->nlambda);
+ snew((dfhist->accum_m)[i],dfhist->nlambda);
+ snew((dfhist->accum_p2)[i],dfhist->nlambda);
+ snew((dfhist->accum_m2)[i],dfhist->nlambda);
+ }
+}
+
+extern void copy_df_history(df_history_t *df_dest, df_history_t *df_source)
+{
+ int i,j;
+
+ init_df_history(df_dest,df_source->nlambda,df_source->wl_delta);
+ df_dest->nlambda = df_source->nlambda;
+ df_dest->bEquil = df_source->bEquil;
+ for (i=0;i<df_dest->nlambda;i++)
+ {
+ df_dest->sum_weights[i] = df_source->sum_weights[i];
+ df_dest->sum_dg[i] = df_source->sum_dg[i];
+ df_dest->sum_minvar[i] = df_source->sum_minvar[i];
+ df_dest->sum_variance[i] = df_source->sum_variance[i];
+ df_dest->n_at_lam[i] = df_source->n_at_lam[i];
+ df_dest->wl_histo[i] = df_source->wl_histo[i];
+ df_dest->accum_p[i] = df_source->accum_p[i];
+ df_dest->accum_m[i] = df_source->accum_m[i];
+ df_dest->accum_p2[i] = df_source->accum_p2[i];
+ df_dest->accum_m2[i] = df_source->accum_m2[i];
+ }
+
+ for (i=0;i<df_dest->nlambda;i++)
+ {
+ for (j=0;j<df_dest->nlambda;j++)
+ {
+ df_dest->Tij[i][j] = df_source->Tij[i][j];
+ df_dest->Tij_empirical[i][j] = df_source->Tij_empirical[i][j];
+ }
+ }
+}
integer less than, equal to, or greater than 0 if \p s1 less than,
identical to, or greater than \p s2.
*/
- static gmx_bool isStringEqNCase(const string s1, const string s2)
+ static gmx_bool isStringEqNCase(const string& s1, const string& s2)
{
return (gmx_strncasecmp(s1.c_str(), s2.c_str(), max(s1.length(), s2.length())) == 0);
}
if (pluginDir != NULL && *pluginDir != '\0')
{
loadedPlugins = Platform::loadPluginsFromDirectory(pluginDir);
- if (loadedPlugins.size() > 0)
+ if (!loadedPlugins.empty())
{
hasLoadedPlugins = true;
usedPluginDir = pluginDir;
if (!hasLoadedPlugins)
{
loadedPlugins = Platform::loadPluginsFromDirectory(OPENMM_PLUGIN_DIR);
- if (loadedPlugins.size() > 0)
+ if (!loadedPlugins.empty())
{
hasLoadedPlugins = true;
usedPluginDir = OPENMM_PLUGIN_DIR;
if (!hasLoadedPlugins)
{
loadedPlugins = Platform::loadPluginsFromDirectory(Platform::getDefaultPluginsDirectory());
- if (loadedPlugins.size() > 0)
+ if (!loadedPlugins.empty())
{
hasLoadedPlugins = true;
usedPluginDir = Platform::getDefaultPluginsDirectory();
int atom3 = ubAtoms[offset++];
/* ubBondForce->addBond(atom1, atom3, */
bondForce->addBond(atom1, atom3,
- idef.iparams[type].u_b.r13, idef.iparams[type].u_b.kUB);
+ idef.iparams[type].u_b.r13A, idef.iparams[type].u_b.kUBA);
/* ubAngleForce->addAngle(atom1, atom2, atom3, */
angleForce->addAngle(atom1, atom2, atom3,
- idef.iparams[type].u_b.theta*M_PI/180.0, idef.iparams[type].u_b.ktheta);
+ idef.iparams[type].u_b.thetaA*M_PI/180.0, idef.iparams[type].u_b.kthetaA);
}
/* Set proper dihedral terms */
return select_res(eargNR,resnr,lh,expl,"ARGININE",nrr,rr);
}
-static const char *get_cystp(int resnr,int nrr,const rtprename_t *rr)
-{
- enum { ecys, ecysH, ecysNR };
- const char *lh[ecysNR] = { "CYS2", "CYS" };
- const char *expl[ecysNR] = {
- "Cysteine in disulfide bridge",
- "Protonated"
- };
-
- return select_res(ecysNR,resnr,lh,expl,"CYSTEINE",nrr,rr);
-
-}
-
static const char *get_histp(int resnr,int nrr,const rtprename_t *rr)
{
const char *expl[ehisNR] = {
for(i=0; i<nrrn; i++) {
fp = fflib_open(rrn[i]);
read_rtprename(rrn[i],fp,&nrtprename,&rtprename);
- fclose(fp);
+ ffclose(fp);
sfree(rrn[i]);
}
sfree(rrn);
}
/* Generate Hydrogen atoms (and termini) in the sequence */
+ printf("Generating any missing hydrogen atoms and/or adding termini.\n");
natom=add_h(&pdba,&x,nah,ah,
cc->nterpairs,cc->ntdb,cc->ctdb,cc->r_start,cc->r_end,bAllowMissing,
NULL,NULL,TRUE,FALSE);
#define MAXPTR 254
#define NOGID 255
+#define MAXLAMBDAS 1024
/* Resource parameters
* Do not change any of these until you read the instruction
energy[STRLEN],user1[STRLEN],user2[STRLEN],vcm[STRLEN],xtc_grps[STRLEN],
couple_moltype[STRLEN],orirefitgrp[STRLEN],egptable[STRLEN],egpexcl[STRLEN],
wall_atomtype[STRLEN],wall_density[STRLEN],deform[STRLEN],QMMM[STRLEN];
-static char foreign_lambda[STRLEN];
+static char fep_lambda[efptNR][STRLEN];
+static char lambda_weights[STRLEN];
static char **pull_grp;
+static char **rot_grp;
static char anneal[STRLEN],anneal_npoints[STRLEN],
anneal_time[STRLEN],anneal_temp[STRLEN];
static char QMmethod[STRLEN],QMbasis[STRLEN],QMcharge[STRLEN],QMmult[STRLEN],
{
snew(opts->include,STRLEN);
snew(opts->define,STRLEN);
+ snew(ir->fepvals,1);
+ snew(ir->expandedvals,1);
+ snew(ir->simtempvals,1);
}
+static void GetSimTemps(int ntemps, t_simtemp *simtemp, double *temperature_lambdas)
+{
+
+ int i;
+
+ for (i=0;i<ntemps;i++)
+ {
+ /* simple linear scaling -- allows more control */
+ if (simtemp->eSimTempScale == esimtempLINEAR)
+ {
+ simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*temperature_lambdas[i];
+ }
+ else if (simtemp->eSimTempScale == esimtempGEOMETRIC) /* should give roughly equal acceptance for constant heat capacity . . . */
+ {
+ simtemp->temperatures[i] = simtemp->simtemp_low * pow(simtemp->simtemp_high/simtemp->simtemp_low,(1.0*i)/(ntemps-1));
+ }
+ else if (simtemp->eSimTempScale == esimtempEXPONENTIAL)
+ {
+ simtemp->temperatures[i] = simtemp->simtemp_low + (simtemp->simtemp_high-simtemp->simtemp_low)*((exp(temperature_lambdas[i])-1)/(exp(1.0)-1));
+ }
+ else
+ {
+ char errorstr[128];
+ sprintf(errorstr,"eSimTempScale=%d not defined",simtemp->eSimTempScale);
+ gmx_fatal(FARGS,errorstr);
+ }
+ }
+}
+
+
+
static void _low_check(gmx_bool b,char *s,warninp_t wi)
{
if (b)
*/
#define CHECK(b) _low_check(b,err_buf,wi)
char err_buf[256],warn_buf[STRLEN];
+ int i,j;
int ns_type=0;
+ real dt_coupl=0;
real dt_pcoupl;
+ int nstcmin;
+ t_lambda *fep = ir->fepvals;
+ t_expanded *expand = ir->expandedvals;
set_warning_line(wi,mdparin,-1);
/* BASIC CUT-OFF STUFF */
+ if (ir->rcoulomb < 0)
+ {
+ warning_error(wi,"rcoulomb should be >= 0");
+ }
+ if (ir->rvdw < 0)
+ {
+ warning_error(wi,"rvdw should be >= 0");
+ }
+ if (ir->rlist < 0)
+ {
+ warning_error(wi,"rlist should be >= 0");
+ }
if (ir->rlist == 0 ||
!((EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > ir->rlist) ||
(EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > ir->rlist))) {
{
ir->etc = etcNO;
}
+ if (ir->eI == eiVVAK) {
+ sprintf(warn_buf,"Integrator method %s is implemented primarily for validation purposes; for molecular dynamics, you should probably be using %s or %s",ei_names[eiVVAK],ei_names[eiMD],ei_names[eiVV]);
+ warning_note(wi,warn_buf);
+ }
if (!EI_DYNAMICS(ir->eI))
{
ir->epc = epcNO;
{
/* nstdhdl should be a multiple of nstcalcenergy */
check_nst("nstcalcenergy",ir->nstcalcenergy,
- "nstdhdl",&ir->nstdhdl,wi);
+ "nstdhdl",&ir->fepvals->nstdhdl,wi);
}
}
}
/* SHAKE / LINCS */
if ( (opts->nshake > 0) && (opts->bMorse) ) {
- sprintf(warn_buf,
- "Using morse bond-potentials while constraining bonds is useless");
- warning(wi,warn_buf);
+ sprintf(warn_buf,
+ "Using morse bond-potentials while constraining bonds is useless");
+ warning(wi,warn_buf);
}
-
- sprintf(err_buf,"shake_tol must be > 0 instead of %g while using shake",
- ir->shake_tol);
- CHECK(((ir->shake_tol <= 0.0) && (opts->nshake>0) &&
- (ir->eConstrAlg == econtSHAKE)));
-
+
+ if ((EI_SD(ir->eI) || ir->eI == eiBD) &&
+ ir->bContinuation && ir->ld_seed != -1) {
+ warning_note(wi,"You are doing a continuation with SD or BD, make sure that ld_seed is different from the previous run (using ld_seed=-1 will ensure this)");
+ }
+ /* verify simulated tempering options */
+
+ if (ir->bSimTemp) {
+ gmx_bool bAllTempZero = TRUE;
+ for (i=0;i<fep->n_lambda;i++)
+ {
+ sprintf(err_buf,"Entry %d for %s must be between 0 and 1, instead is %g",i,efpt_names[efptTEMPERATURE],fep->all_lambda[efptTEMPERATURE][i]);
+ CHECK((fep->all_lambda[efptTEMPERATURE][i] < 0) || (fep->all_lambda[efptTEMPERATURE][i] > 1));
+ if (fep->all_lambda[efptTEMPERATURE][i] > 0)
+ {
+ bAllTempZero = FALSE;
+ }
+ }
+ sprintf(err_buf,"if simulated tempering is on, temperature-lambdas may not be all zero");
+ CHECK(bAllTempZero==TRUE);
+
+ sprintf(err_buf,"Simulated tempering is currently only compatible with md-vv");
+ CHECK(ir->eI != eiVV);
+
+ /* check compatability of the temperature coupling with simulated tempering */
+
+ if (ir->etc == etcNOSEHOOVER) {
+ sprintf(warn_buf,"Nose-Hoover based temperature control such as [%s] my not be entirelyconsistent with simulated tempering",etcoupl_names[ir->etc]);
+ warning_note(wi,warn_buf);
+ }
+
+ /* check that the temperatures make sense */
+
+ sprintf(err_buf,"Higher simulated tempering temperature (%g) must be >= than the simulated tempering lower temperature (%g)",ir->simtempvals->simtemp_high,ir->simtempvals->simtemp_low);
+ CHECK(ir->simtempvals->simtemp_high <= ir->simtempvals->simtemp_low);
+
+ sprintf(err_buf,"Higher simulated tempering temperature (%g) must be >= zero",ir->simtempvals->simtemp_high);
+ CHECK(ir->simtempvals->simtemp_high <= 0);
+
+ sprintf(err_buf,"Lower simulated tempering temperature (%g) must be >= zero",ir->simtempvals->simtemp_low);
+ CHECK(ir->simtempvals->simtemp_low <= 0);
+ }
+
+ /* verify free energy options */
+
+ if (ir->efep != efepNO) {
+ fep = ir->fepvals;
+ sprintf(err_buf,"The soft-core power is %d and can only be 1 or 2",
+ fep->sc_power);
+ CHECK(fep->sc_alpha!=0 && fep->sc_power!=1 && fep->sc_power!=2);
+
+ sprintf(err_buf,"The soft-core sc-r-power is %d and can only be 6 or 48",
+ (int)fep->sc_r_power);
+ CHECK(fep->sc_alpha!=0 && fep->sc_r_power!=6.0 && fep->sc_r_power!=48.0);
+
+ /* check validity of options */
+ if (fep->n_lambda > 0 && ir->rlist < max(ir->rvdw,ir->rcoulomb))
+ {
+ sprintf(warn_buf,
+ "For foreign lambda free energy differences it is assumed that the soft-core interactions have no effect beyond the neighborlist cut-off");
+ warning(wi,warn_buf);
+ }
+
+ sprintf(err_buf,"Can't use postive delta-lambda (%g) if initial state/lambda does not start at zero",fep->delta_lambda);
+ CHECK(fep->delta_lambda > 0 && ((fep->init_fep_state !=0) || (fep->init_lambda !=0)));
+
+ sprintf(err_buf,"Can't use postive delta-lambda (%g) with expanded ensemble simulations",fep->delta_lambda);
+ CHECK(fep->delta_lambda > 0 && (ir->efep == efepEXPANDED));
+
+ sprintf(err_buf,"Free-energy not implemented for Ewald");
+ CHECK(ir->coulombtype==eelEWALD);
+
+ /* check validty of lambda inputs */
+ sprintf(err_buf,"initial thermodynamic state %d does not exist, only goes to %d",fep->init_fep_state,fep->n_lambda);
+ CHECK((fep->init_fep_state > fep->n_lambda));
+
+ for (j=0;j<efptNR;j++)
+ {
+ for (i=0;i<fep->n_lambda;i++)
+ {
+ sprintf(err_buf,"Entry %d for %s must be between 0 and 1, instead is %g",i,efpt_names[j],fep->all_lambda[j][i]);
+ CHECK((fep->all_lambda[j][i] < 0) || (fep->all_lambda[j][i] > 1));
+ }
+ }
+
+ if ((fep->sc_alpha>0) && (!fep->bScCoul))
+ {
+ for (i=0;i<fep->n_lambda;i++)
+ {
+ sprintf(err_buf,"For state %d, vdw-lambdas (%f) is changing with vdw softcore, while coul-lambdas (%f) is nonzero without coulomb softcore: this will lead to crashes, and is not supported.",i,fep->all_lambda[efptVDW][i],
+ fep->all_lambda[efptCOUL][i]);
+ CHECK((fep->sc_alpha>0) &&
+ (((fep->all_lambda[efptCOUL][i] > 0.0) &&
+ (fep->all_lambda[efptCOUL][i] < 1.0)) &&
+ ((fep->all_lambda[efptVDW][i] > 0.0) &&
+ (fep->all_lambda[efptVDW][i] < 1.0))));
+ }
+ }
+
+ if ((fep->bScCoul) && (EEL_PME(ir->coulombtype)))
+ {
+ sprintf(warn_buf,"With coulomb soft core, the reciprocal space calculation will not necessarily cancel. It may be necessary to decrease the reciprocal space energy, and increase the cutoff radius to get sufficiently close matches to energies with free energy turned off.");
+ warning(wi, warn_buf);
+ }
+
+ /* Free Energy Checks -- In an ideal world, slow growth and FEP would
+ be treated differently, but that's the next step */
+
+ for (i=0;i<efptNR;i++) {
+ for (j=0;j<fep->n_lambda;j++) {
+ sprintf(err_buf,"%s[%d] must be between 0 and 1",efpt_names[i],j);
+ CHECK((fep->all_lambda[i][j] < 0) || (fep->all_lambda[i][j] > 1));
+ }
+ }
+ }
+
+ if ((ir->bSimTemp) || (ir->efep == efepEXPANDED)) {
+ fep = ir->fepvals;
+ expand = ir->expandedvals;
+
+ /* checking equilibration of weights inputs for validity */
+
+ sprintf(err_buf,"weight-equil-number-all-lambda (%d) is ignored if lmc-weights-equil is not equal to %s",
+ expand->equil_n_at_lam,elmceq_names[elmceqNUMATLAM]);
+ CHECK((expand->equil_n_at_lam>0) && (expand->elmceq!=elmceqNUMATLAM));
+
+ sprintf(err_buf,"weight-equil-number-samples (%d) is ignored if lmc-weights-equil is not equal to %s",
+ expand->equil_samples,elmceq_names[elmceqSAMPLES]);
+ CHECK((expand->equil_samples>0) && (expand->elmceq!=elmceqSAMPLES));
+
+ sprintf(err_buf,"weight-equil-number-steps (%d) is ignored if lmc-weights-equil is not equal to %s",
+ expand->equil_steps,elmceq_names[elmceqSTEPS]);
+ CHECK((expand->equil_steps>0) && (expand->elmceq!=elmceqSTEPS));
+
+ sprintf(err_buf,"weight-equil-wl-delta (%d) is ignored if lmc-weights-equil is not equal to %s",
+ expand->equil_samples,elmceq_names[elmceqWLDELTA]);
+ CHECK((expand->equil_wl_delta>0) && (expand->elmceq!=elmceqWLDELTA));
+
+ sprintf(err_buf,"weight-equil-count-ratio (%f) is ignored if lmc-weights-equil is not equal to %s",
+ expand->equil_ratio,elmceq_names[elmceqRATIO]);
+ CHECK((expand->equil_ratio>0) && (expand->elmceq!=elmceqRATIO));
+
+ sprintf(err_buf,"weight-equil-number-all-lambda (%d) must be a positive integer if lmc-weights-equil=%s",
+ expand->equil_n_at_lam,elmceq_names[elmceqNUMATLAM]);
+ CHECK((expand->equil_n_at_lam<=0) && (expand->elmceq==elmceqNUMATLAM));
+
+ sprintf(err_buf,"weight-equil-number-samples (%d) must be a positive integer if lmc-weights-equil=%s",
+ expand->equil_samples,elmceq_names[elmceqSAMPLES]);
+ CHECK((expand->equil_samples<=0) && (expand->elmceq==elmceqSAMPLES));
+
+ sprintf(err_buf,"weight-equil-number-steps (%d) must be a positive integer if lmc-weights-equil=%s",
+ expand->equil_steps,elmceq_names[elmceqSTEPS]);
+ CHECK((expand->equil_steps<=0) && (expand->elmceq==elmceqSTEPS));
+
+ sprintf(err_buf,"weight-equil-wl-delta (%f) must be > 0 if lmc-weights-equil=%s",
+ expand->equil_wl_delta,elmceq_names[elmceqWLDELTA]);
+ CHECK((expand->equil_wl_delta<=0) && (expand->elmceq==elmceqWLDELTA));
+
+ sprintf(err_buf,"weight-equil-count-ratio (%f) must be > 0 if lmc-weights-equil=%s",
+ expand->equil_ratio,elmceq_names[elmceqRATIO]);
+ CHECK((expand->equil_ratio<=0) && (expand->elmceq==elmceqRATIO));
+
+ sprintf(err_buf,"lmc-weights-equil=%s only possible when lmc-stats = %s or lmc-stats %s",
+ elmceq_names[elmceqWLDELTA],elamstats_names[elamstatsWL],elamstats_names[elamstatsWWL]);
+ CHECK((expand->elmceq==elmceqWLDELTA) && (!EWL(expand->elamstats)));
+
+ sprintf(err_buf,"lmc-repeats (%d) must be greater than 0",expand->lmc_repeats);
+ CHECK((expand->lmc_repeats <= 0));
+ sprintf(err_buf,"minimum-var-min (%d) must be greater than 0",expand->minvarmin);
+ CHECK((expand->minvarmin <= 0));
+ sprintf(err_buf,"weight-c-range (%d) must be greater or equal to 0",expand->c_range);
+ CHECK((expand->c_range < 0));
+ sprintf(err_buf,"init-lambda-state (%d) must be zero if lmc-forced-nstart (%d)> 0 and lmc-move != 'no'",
+ fep->init_fep_state, expand->lmc_forced_nstart);
+ CHECK((fep->init_fep_state!=0) && (expand->lmc_forced_nstart>0) && (expand->elmcmove!=elmcmoveNO));
+ sprintf(err_buf,"lmc-forced-nstart (%d) must not be negative",expand->lmc_forced_nstart);
+ CHECK((expand->lmc_forced_nstart < 0));
+ sprintf(err_buf,"init-lambda-state (%d) must be in the interval [0,number of lambdas)",fep->init_fep_state);
+ CHECK((fep->init_fep_state < 0) || (fep->init_fep_state >= fep->n_lambda));
+
+ sprintf(err_buf,"init-wl-delta (%f) must be greater than or equal to 0",expand->init_wl_delta);
+ CHECK((expand->init_wl_delta < 0));
+ sprintf(err_buf,"wl-ratio (%f) must be between 0 and 1",expand->wl_ratio);
+ CHECK((expand->wl_ratio <= 0) || (expand->wl_ratio >= 1));
+ sprintf(err_buf,"wl-scale (%f) must be between 0 and 1",expand->wl_scale);
+ CHECK((expand->wl_scale <= 0) || (expand->wl_scale >= 1));
+
+ /* if there is no temperature control, we need to specify an MC temperature */
+ sprintf(err_buf,"If there is no temperature control, and lmc-mcmove!= 'no',mc_temperature must be set to a positive number");
+ if (expand->nstTij > 0)
+ {
+ sprintf(err_buf,"nst-transition-matrix (%d) must be an integer multiple of nstlog (%d)",
+ expand->nstTij,ir->nstlog);
+ CHECK((mod(expand->nstTij,ir->nstlog)!=0));
+ }
+ }
+
/* PBC/WALLS */
sprintf(err_buf,"walls only work with pbc=%s",epbc_names[epbcXY]);
CHECK(ir->nwall && ir->ePBC!=epbcXY);
}
sprintf(err_buf,"Can not have Ewald with pbc=%s",epbc_names[ir->ePBC]);
CHECK(EEL_FULL(ir->coulombtype));
-
+
sprintf(err_buf,"Can not have dispersion correction with pbc=%s",
epbc_names[ir->ePBC]);
CHECK(ir->eDispCorr != edispcNO);
"rcoulomb and rvdw set to zero",
eel_names[eelCUT],eel_names[eelUSER],epbc_names[epbcNONE]);
CHECK(((ir->coulombtype != eelCUT) && (ir->coulombtype != eelUSER)) ||
- (ir->ePBC != epbcNONE) ||
+ (ir->ePBC != epbcNONE) ||
(ir->rcoulomb != 0.0) || (ir->rvdw != 0.0));
if (ir->nstlist < 0) {
warning(wi,"If you want to remove the rotation around the center of mass, you should set comm_mode = Angular instead of setting nstcomm < 0. nstcomm is modified to its absolute value");
ir->nstcomm = abs(ir->nstcomm);
}
-
+
if (ir->nstcalcenergy > 0 && ir->nstcomm < ir->nstcalcenergy) {
warning_note(wi,"nstcomm < nstcalcenergy defeats the purpose of nstcalcenergy, setting nstcomm to nstcalcenergy");
- ir->nstcomm = ir->nstcalcenergy;
+ ir->nstcomm = ir->nstcalcenergy;
}
if (ir->comm_mode == ecmANGULAR) {
warning(wi,"Removing the rotation around the center of mass in a periodic system (this is not a problem when you have only one molecule).");
}
}
-
+
if (EI_STATE_VELOCITY(ir->eI) && ir->ePBC == epbcNONE && ir->comm_mode != ecmANGULAR) {
warning_note(wi,"Tumbling and or flying ice-cubes: We are not removing rotation around center of mass in a non-periodic system. You should probably set comm_mode = ANGULAR.");
}
- sprintf(err_buf,"Free-energy not implemented for Ewald and PPPM");
- CHECK((ir->coulombtype==eelEWALD || ir->coulombtype==eelPPPM)
- && (ir->efep!=efepNO));
-
sprintf(err_buf,"Twin-range neighbour searching (NS) with simple NS"
" algorithm not implemented");
- CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
+ CHECK(((ir->rcoulomb > ir->rlist) || (ir->rvdw > ir->rlist))
&& (ir->ns_type == ensSIMPLE));
-
- /* TEMPERATURE COUPLING */
- if (ir->etc == etcYES)
+
+ /* TEMPERATURE COUPLING */
+ if (ir->etc == etcYES)
{
ir->etc = etcBERENDSEN;
warning_note(wi,"Old option for temperature coupling given: "
"changing \"yes\" to \"Berendsen\"\n");
}
-
- if (ir->etc == etcNOSEHOOVER)
+
+ if ((ir->etc == etcNOSEHOOVER) || (ir->epc == epcMTTK))
{
- if (ir->opts.nhchainlength < 1)
+ if (ir->opts.nhchainlength < 1)
{
sprintf(warn_buf,"number of Nose-Hoover chains (currently %d) cannot be less than 1,reset to 1\n",ir->opts.nhchainlength);
ir->opts.nhchainlength =1;
ir->opts.nhchainlength = 0;
}
+ if (ir->eI == eiVVAK) {
+ sprintf(err_buf,"%s implemented primarily for validation, and requires nsttcouple = 1 and nstpcouple = 1.",
+ ei_names[eiVVAK]);
+ CHECK((ir->nsttcouple != 1) || (ir->nstpcouple != 1));
+ }
+
+ if (ETC_ANDERSEN(ir->etc))
+ {
+ sprintf(err_buf,"%s temperature control not supported for integrator %s.",etcoupl_names[ir->etc],ei_names[ir->eI]);
+ CHECK(!(EI_VV(ir->eI)));
+
+ for (i=0;i<ir->opts.ngtc;i++)
+ {
+ sprintf(err_buf,"all tau_t must currently be equal using Andersen temperature control, violated for group %d",i);
+ CHECK(ir->opts.tau_t[0] != ir->opts.tau_t[i]);
+ sprintf(err_buf,"all tau_t must be postive using Andersen temperature control, tau_t[%d]=%10.6f",
+ i,ir->opts.tau_t[i]);
+ CHECK(ir->opts.tau_t[i]<0);
+ }
+ if (ir->nstcomm > 0 && (ir->etc == etcANDERSEN)) {
+ sprintf(warn_buf,"Center of mass removal not necessary for %s. All velocities of coupled groups are rerandomized periodically, so flying ice cube errors will not occur.",etcoupl_names[ir->etc]);
+ warning_note(wi,warn_buf);
+ }
+
+ sprintf(err_buf,"nstcomm must be 1, not %d for %s, as velocities of atoms in coupled groups are randomized every time step",ir->nstcomm,etcoupl_names[ir->etc]);
+ CHECK(ir->nstcomm > 1 && (ir->etc == etcANDERSEN));
+
+ for (i=0;i<ir->opts.ngtc;i++)
+ {
+ int nsteps = (int)(ir->opts.tau_t[i]/ir->delta_t);
+ sprintf(err_buf,"tau_t/delta_t for group %d for temperature control method %s must be a multiple of nstcomm (%d), as velocities of atoms in coupled groups are randomized every time step. The input tau_t (%8.3f) leads to %d steps per randomization",i,etcoupl_names[ir->etc],ir->nstcomm,ir->opts.tau_t[i],nsteps);
+ CHECK((nsteps % ir->nstcomm) && (ir->etc == etcANDERSENMASSIVE));
+ }
+ }
if (ir->etc == etcBERENDSEN)
{
sprintf(warn_buf,"The %s thermostat does not generate the correct kinetic energy distribution. You might want to consider using the %s thermostat.",
warning_note(wi,warn_buf);
}
- if ((ir->etc==etcNOSEHOOVER || ir->etc==etcANDERSEN || ir->etc==etcANDERSENINTERVAL)
+ if ((ir->etc==etcNOSEHOOVER || ETC_ANDERSEN(ir->etc))
&& ir->epc==epcBERENDSEN)
{
sprintf(warn_buf,"Using Berendsen pressure coupling invalidates the "
{
dt_pcoupl = ir->nstpcouple*ir->delta_t;
- sprintf(err_buf,"tau_p must be > 0 instead of %g\n",ir->tau_p);
+ sprintf(err_buf,"tau-p must be > 0 instead of %g\n",ir->tau_p);
CHECK(ir->tau_p <= 0);
-
+
if (ir->tau_p/dt_pcoupl < pcouple_min_integration_steps(ir->epc))
{
- sprintf(warn_buf,"For proper integration of the %s barostat, tau_p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
+ sprintf(warn_buf,"For proper integration of the %s barostat, tau-p (%g) should be at least %d times larger than nstpcouple*dt (%g)",
EPCOUPLTYPE(ir->epc),ir->tau_p,pcouple_min_integration_steps(ir->epc),dt_pcoupl);
warning(wi,warn_buf);
- }
-
- sprintf(err_buf,"compressibility must be > 0 when using pressure"
+ }
+
+ sprintf(err_buf,"compressibility must be > 0 when using pressure"
" coupling %s\n",EPCOUPLTYPE(ir->epc));
- CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
- ir->compress[ZZ][ZZ] < 0 ||
+ CHECK(ir->compress[XX][XX] < 0 || ir->compress[YY][YY] < 0 ||
+ ir->compress[ZZ][ZZ] < 0 ||
(trace(ir->compress) == 0 && ir->compress[YY][XX] <= 0 &&
ir->compress[ZZ][XX] <= 0 && ir->compress[ZZ][YY] <= 0));
- if (epcPARRINELLORAHMAN == ir->epct && opts->bGenVel)
- sprintf(err_buf,"pressure coupling with PPPM not implemented, use PME");
- CHECK(ir->coulombtype == eelPPPM);
-
+ if (epcPARRINELLORAHMAN == ir->epc && opts->bGenVel)
{
sprintf(warn_buf,
"You are generating velocities so I am assuming you "
"are equilibrating a system. You are using "
- "Parrinello-Rahman pressure coupling, but this can be "
+ "%s pressure coupling, but this can be "
"unstable for equilibration. If your system crashes, try "
"equilibrating first with Berendsen pressure coupling. If "
"you are not equilibrating the system, you can probably "
- "ignore this warning.");
+ "ignore this warning.",
+ epcoupl_names[ir->epc]);
warning(wi,warn_buf);
}
}
- else if (ir->coulombtype == eelPPPM)
- {
- sprintf(warn_buf,"The pressure with PPPM is incorrect, if you need the pressure use PME");
- warning(wi,warn_buf);
- }
-
+
if (EI_VV(ir->eI))
{
if (ir->epc > epcNO)
{
- if (ir->epc!=epcMTTK)
+ if ((ir->epc!=epcBERENDSEN) && (ir->epc!=epcMTTK))
{
- warning_error(wi,"NPT only defined for vv using Martyna-Tuckerman-Tobias-Klein equations");
+ warning_error(wi,"for md-vv and md-vv-avek, can only use Berendsen and Martyna-Tuckerman-Tobias-Klein (MTTK) equations for pressure control; MTTK is equivalent to Parrinello-Rahman.");
}
}
}
/* ELECTROSTATICS */
/* More checks are in triple check (grompp.c) */
- if (ir->coulombtype == eelPPPM)
- {
- warning_error(wi,"PPPM is not functional in the current version, we plan to implement PPPM through a small modification of the PME code");
- }
if (ir->coulombtype == eelSWITCH) {
sprintf(warn_buf,"coulombtype = %s is only for testing purposes and can lead to serious artifacts, advice: use coulombtype = %s",
}
if (ir->epsilon_r!=1 && ir->implicit_solvent==eisGBSA) {
- sprintf(warn_buf,"epsilon_r = %g with GB implicit solvent, will use this value for inner dielectric",ir->epsilon_r);
+ sprintf(warn_buf,"epsilon-r = %g with GB implicit solvent, will use this value for inner dielectric",ir->epsilon_r);
warning_note(wi,warn_buf);
}
if (EEL_RF(ir->coulombtype) && ir->epsilon_rf==1 && ir->epsilon_r!=1) {
- sprintf(warn_buf,"epsilon_r = %g and epsilon_rf = 1 with reaction field, assuming old format and exchanging epsilon_r and epsilon_rf",ir->epsilon_r);
+ sprintf(warn_buf,"epsilon-r = %g and epsilon-rf = 1 with reaction field, assuming old format and exchanging epsilon-r and epsilon-rf",ir->epsilon_r);
warning(wi,warn_buf);
ir->epsilon_rf = ir->epsilon_r;
ir->epsilon_r = 1.0;
}
if (getenv("GALACTIC_DYNAMICS") == NULL) {
- sprintf(err_buf,"epsilon_r must be >= 0 instead of %g\n",ir->epsilon_r);
+ sprintf(err_buf,"epsilon-r must be >= 0 instead of %g\n",ir->epsilon_r);
CHECK(ir->epsilon_r < 0);
}
/* reaction field (at the cut-off) */
if (ir->coulombtype == eelRF_ZERO) {
- sprintf(err_buf,"With coulombtype = %s, epsilon_rf must be 0",
+ sprintf(err_buf,"With coulombtype = %s, epsilon-rf must be 0",
eel_names[ir->coulombtype]);
CHECK(ir->epsilon_rf != 0);
}
- sprintf(err_buf,"epsilon_rf must be >= epsilon_r");
+ sprintf(err_buf,"epsilon-rf must be >= epsilon-r");
CHECK((ir->epsilon_rf < ir->epsilon_r && ir->epsilon_rf != 0) ||
(ir->epsilon_r == 0));
if (ir->epsilon_rf == ir->epsilon_r) {
- sprintf(warn_buf,"Using epsilon_rf = epsilon_r with %s does not make sense",
+ sprintf(warn_buf,"Using epsilon-rf = epsilon-r with %s does not make sense",
eel_names[ir->coulombtype]);
warning(wi,warn_buf);
}
eel_names[ir->coulombtype]);
CHECK(ir->rcoulomb > ir->rlist);
} else {
- if (ir->coulombtype == eelPME) {
+ if (ir->coulombtype == eelPME || ir->coulombtype == eelP3M_AD) {
sprintf(err_buf,
"With coulombtype = %s, rcoulomb must be equal to rlist\n"
"If you want optimal energy conservation or exact integration use %s",
if (EEL_PME(ir->coulombtype)) {
if (ir->pme_order < 3) {
- warning_error(wi,"pme_order can not be smaller than 3");
+ warning_error(wi,"pme-order can not be smaller than 3");
}
}
if (ir->nwall==2 && EEL_FULL(ir->coulombtype)) {
if (ir->ewald_geometry == eewg3D) {
- sprintf(warn_buf,"With pbc=%s you should use ewald_geometry=%s",
+ sprintf(warn_buf,"With pbc=%s you should use ewald-geometry=%s",
epbc_names[ir->ePBC],eewg_names[eewg3DC]);
warning(wi,warn_buf);
}
/* This check avoids extra pbc coding for exclusion corrections */
- sprintf(err_buf,"wall_ewald_zfac should be >= 2");
+ sprintf(err_buf,"wall-ewald-zfac should be >= 2");
CHECK(ir->wall_ewald_zfac < 2);
}
if (EVDW_SWITCHED(ir->vdwtype)) {
- sprintf(err_buf,"With vdwtype = %s rvdw_switch must be < rvdw",
+ sprintf(err_buf,"With vdwtype = %s rvdw-switch must be < rvdw",
evdw_names[ir->vdwtype]);
CHECK(ir->rvdw_switch >= ir->rvdw);
} else if (ir->vdwtype == evdwCUT) {
warning(wi,"Using L-BFGS with nbfgscorr<=0 just gets you steepest descent.");
}
- /* FREE ENERGY */
- if (ir->efep != efepNO) {
- sprintf(err_buf,"The soft-core power is %d and can only be 1 or 2",
- ir->sc_power);
- CHECK(ir->sc_alpha!=0 && ir->sc_power!=1 && ir->sc_power!=2);
+ /* ENERGY CONSERVATION */
+ if (ir_NVE(ir))
+ {
+ if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0)
+ {
+ sprintf(warn_buf,"You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
+ evdw_names[evdwSHIFT]);
+ warning_note(wi,warn_buf);
+ }
+ if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0)
+ {
+ sprintf(warn_buf,"You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
+ eel_names[eelPMESWITCH],eel_names[eelRF_ZERO]);
+ warning_note(wi,warn_buf);
+ }
}
- /* ENERGY CONSERVATION */
- if (ir_NVE(ir))
- {
- if (!EVDW_MIGHT_BE_ZERO_AT_CUTOFF(ir->vdwtype) && ir->rvdw > 0)
- {
- sprintf(warn_buf,"You are using a cut-off for VdW interactions with NVE, for good energy conservation use vdwtype = %s (possibly with DispCorr)",
- evdw_names[evdwSHIFT]);
- warning_note(wi,warn_buf);
- }
- if (!EEL_MIGHT_BE_ZERO_AT_CUTOFF(ir->coulombtype) && ir->rcoulomb > 0)
- {
- sprintf(warn_buf,"You are using a cut-off for electrostatics with NVE, for good energy conservation use coulombtype = %s or %s",
- eel_names[eelPMESWITCH],eel_names[eelRF_ZERO]);
- warning_note(wi,warn_buf);
- }
- }
-
/* IMPLICIT SOLVENT */
if(ir->coulombtype==eelGB_NOTUSED)
{
ir->implicit_solvent=eisGBSA;
fprintf(stderr,"Note: Old option for generalized born electrostatics given:\n"
"Changing coulombtype from \"generalized-born\" to \"cut-off\" and instead\n"
- "setting implicit_solvent value to \"GBSA\" in input section.\n");
+ "setting implicit-solvent value to \"GBSA\" in input section.\n");
}
if(ir->sa_algorithm==esaSTILL)
}
}
+
+ if (ir->bAdress && !EI_SD(ir->eI)){
+ warning_error(wi,"AdresS simulation supports only stochastic dynamics");
+ }
+ if (ir->bAdress && ir->epc != epcNO){
+ warning_error(wi,"AdresS simulation does not support pressure coupling");
+ }
+ if (ir->bAdress && (EEL_FULL(ir->coulombtype))){
+ warning_error(wi,"AdresS simulation does not support long-range electrostatics");
+ }
+
}
-static int str_nelem(const char *str,int maxptr,char *ptr[])
+/* count the number of text elemets separated by whitespace in a string.
+ str = the input string
+ maxptr = the maximum number of allowed elements
+ ptr = the output array of pointers to the first character of each element
+ returns: the number of elements. */
+int str_nelem(const char *str,int maxptr,char *ptr[])
{
int np=0;
char *copy0,*copy;
return np;
}
-static void parse_n_double(char *str,int *n,double **r)
+/* interpret a number of doubles from a string and put them in an array,
+ after allocating space for them.
+ str = the input string
+ n = the (pre-allocated) number of doubles read
+ r = the output array of doubles. */
+static void parse_n_real(char *str,int *n,real **r)
{
char *ptr[MAXPTR];
int i;
}
}
+static void do_fep_params(t_inputrec *ir, char fep_lambda[][STRLEN],char weights[STRLEN]) {
+
+ int i,j,max_n_lambda,nweights,nfep[efptNR];
+ t_lambda *fep = ir->fepvals;
+ t_expanded *expand = ir->expandedvals;
+ real **count_fep_lambdas;
+ gmx_bool bOneLambda = TRUE;
+
+ snew(count_fep_lambdas,efptNR);
+
+ /* FEP input processing */
+ /* first, identify the number of lambda values for each type.
+ All that are nonzero must have the same number */
+
+ for (i=0;i<efptNR;i++)
+ {
+ parse_n_real(fep_lambda[i],&(nfep[i]),&(count_fep_lambdas[i]));
+ }
+
+ /* now, determine the number of components. All must be either zero, or equal. */
+
+ max_n_lambda = 0;
+ for (i=0;i<efptNR;i++)
+ {
+ if (nfep[i] > max_n_lambda) {
+ max_n_lambda = nfep[i]; /* here's a nonzero one. All of them
+ must have the same number if its not zero.*/
+ break;
+ }
+ }
+
+ for (i=0;i<efptNR;i++)
+ {
+ if (nfep[i] == 0)
+ {
+ ir->fepvals->separate_dvdl[i] = FALSE;
+ }
+ else if (nfep[i] == max_n_lambda)
+ {
+ if (i!=efptTEMPERATURE) /* we treat this differently -- not really a reason to compute the derivative with
+ respect to the temperature currently */
+ {
+ ir->fepvals->separate_dvdl[i] = TRUE;
+ }
+ }
+ else
+ {
+ gmx_fatal(FARGS,"Number of lambdas (%d) for FEP type %s not equal to number of other types (%d)",
+ nfep[i],efpt_names[i],max_n_lambda);
+ }
+ }
+ /* we don't print out dhdl if the temperature is changing, since we can't correctly define dhdl in this case */
+ ir->fepvals->separate_dvdl[efptTEMPERATURE] = FALSE;
+
+ /* the number of lambdas is the number we've read in, which is either zero
+ or the same for all */
+ fep->n_lambda = max_n_lambda;
+
+ /* allocate space for the array of lambda values */
+ snew(fep->all_lambda,efptNR);
+ /* if init_lambda is defined, we need to set lambda */
+ if ((fep->init_lambda > 0) && (fep->n_lambda == 0))
+ {
+ ir->fepvals->separate_dvdl[efptFEP] = TRUE;
+ }
+ /* otherwise allocate the space for all of the lambdas, and transfer the data */
+ for (i=0;i<efptNR;i++)
+ {
+ snew(fep->all_lambda[i],fep->n_lambda);
+ if (nfep[i] > 0) /* if it's zero, then the count_fep_lambda arrays
+ are zero */
+ {
+ for (j=0;j<fep->n_lambda;j++)
+ {
+ fep->all_lambda[i][j] = (double)count_fep_lambdas[i][j];
+ }
+ sfree(count_fep_lambdas[i]);
+ }
+ }
+ sfree(count_fep_lambdas);
+
+ /* "fep-vals" is either zero or the full number. If zero, we'll need to define fep-lambdas for internal
+ bookkeeping -- for now, init_lambda */
+
+ if ((nfep[efptFEP] == 0) && (fep->init_lambda >= 0) && (fep->init_lambda <= 1))
+ {
+ for (i=0;i<fep->n_lambda;i++)
+ {
+ fep->all_lambda[efptFEP][i] = fep->init_lambda;
+ }
+ }
+
+ /* check to see if only a single component lambda is defined, and soft core is defined.
+ In this case, turn on coulomb soft core */
+
+ if (max_n_lambda == 0)
+ {
+ bOneLambda = TRUE;
+ }
+ else
+ {
+ for (i=0;i<efptNR;i++)
+ {
+ if ((nfep[i] != 0) && (i!=efptFEP))
+ {
+ bOneLambda = FALSE;
+ }
+ }
+ }
+ if ((bOneLambda) && (fep->sc_alpha > 0))
+ {
+ fep->bScCoul = TRUE;
+ }
+
+ /* Fill in the others with the efptFEP if they are not explicitly
+ specified (i.e. nfep[i] == 0). This means if fep is not defined,
+ they are all zero. */
+
+ for (i=0;i<efptNR;i++)
+ {
+ if ((nfep[i] == 0) && (i!=efptFEP))
+ {
+ for (j=0;j<fep->n_lambda;j++)
+ {
+ fep->all_lambda[i][j] = fep->all_lambda[efptFEP][j];
+ }
+ }
+ }
+
+
+ /* make it easier if sc_r_power = 48 by increasing it to the 4th power, to be in the right scale. */
+ if (fep->sc_r_power == 48)
+ {
+ if (fep->sc_alpha > 0.1)
+ {
+ gmx_fatal(FARGS,"sc_alpha (%f) for sc_r_power = 48 should usually be between 0.001 and 0.004", fep->sc_alpha);
+ }
+ }
+
+ expand = ir->expandedvals;
+ /* now read in the weights */
+ parse_n_real(weights,&nweights,&(expand->init_lambda_weights));
+ if (nweights == 0)
+ {
+ expand->bInit_weights = FALSE;
+ snew(expand->init_lambda_weights,fep->n_lambda); /* initialize to zero */
+ }
+ else if (nweights != fep->n_lambda)
+ {
+ gmx_fatal(FARGS,"Number of weights (%d) is not equal to number of lambda values (%d)",
+ nweights,fep->n_lambda);
+ }
+ else
+ {
+ expand->bInit_weights = TRUE;
+ }
+ if ((expand->nstexpanded < 0) && (ir->efep != efepNO)) {
+ expand->nstexpanded = fep->nstdhdl;
+ /* if you don't specify nstexpanded when doing expanded ensemble free energy calcs, it is set to nstdhdl */
+ }
+ if ((expand->nstexpanded < 0) && ir->bSimTemp) {
+ expand->nstexpanded = ir->nstlist;
+ /* if you don't specify nstexpanded when doing expanded ensemble simulated tempering, it is set to nstlist*/
+ }
+}
+
+
+static void do_simtemp_params(t_inputrec *ir) {
+
+ snew(ir->simtempvals->temperatures,ir->fepvals->n_lambda);
+ GetSimTemps(ir->fepvals->n_lambda,ir->simtempvals,ir->fepvals->all_lambda[efptTEMPERATURE]);
+
+ return;
+}
+
static void do_wall_params(t_inputrec *ir,
char *wall_atomtype, char *wall_density,
t_gromppopts *opts)
nstr = str_nelem(wall_density,MAXPTR,names);
if (nstr != ir->nwall)
{
- gmx_fatal(FARGS,"Expected %d elements for wall_density, found %d",ir->nwall,nstr);
+ gmx_fatal(FARGS,"Expected %d elements for wall-density, found %d",ir->nwall,nstr);
}
for(i=0; i<ir->nwall; i++)
{
sscanf(names[i],"%lf",&dbl);
if (dbl <= 0)
{
- gmx_fatal(FARGS,"wall_density[%d] = %f\n",i,dbl);
+ gmx_fatal(FARGS,"wall-density[%d] = %f\n",i,dbl);
}
ir->wall_density[i] = dbl;
}
}
}
+void read_expandedparams(int *ninp_p,t_inpfile **inp_p,
+ t_expanded *expand,warninp_t wi)
+{
+ int ninp,nerror=0;
+ t_inpfile *inp;
+
+ ninp = *ninp_p;
+ inp = *inp_p;
+
+ /* read expanded ensemble parameters */
+ CCTYPE ("expanded ensemble variables");
+ ITYPE ("nstexpanded",expand->nstexpanded,-1);
+ EETYPE("lmc-stats", expand->elamstats, elamstats_names);
+ EETYPE("lmc-move", expand->elmcmove, elmcmove_names);
+ EETYPE("lmc-weights-equil",expand->elmceq,elmceq_names);
+ ITYPE ("weight-equil-number-all-lambda",expand->equil_n_at_lam,-1);
+ ITYPE ("weight-equil-number-samples",expand->equil_samples,-1);
+ ITYPE ("weight-equil-number-steps",expand->equil_steps,-1);
+ RTYPE ("weight-equil-wl-delta",expand->equil_wl_delta,-1);
+ RTYPE ("weight-equil-count-ratio",expand->equil_ratio,-1);
+ CCTYPE("Seed for Monte Carlo in lambda space");
+ ITYPE ("lmc-seed",expand->lmc_seed,-1);
+ RTYPE ("mc-temperature",expand->mc_temp,-1);
+ ITYPE ("lmc-repeats",expand->lmc_repeats,1);
+ ITYPE ("lmc-gibbsdelta",expand->gibbsdeltalam,-1);
+ ITYPE ("lmc-forced-nstart",expand->lmc_forced_nstart,0);
+ EETYPE("symmetrized-transition-matrix", expand->bSymmetrizedTMatrix, yesno_names);
+ ITYPE("nst-transition-matrix", expand->nstTij, -1);
+ ITYPE ("mininum-var-min",expand->minvarmin, 100); /*default is reasonable */
+ ITYPE ("weight-c-range",expand->c_range, 0); /* default is just C=0 */
+ RTYPE ("wl-scale",expand->wl_scale,0.8);
+ RTYPE ("wl-ratio",expand->wl_ratio,0.8);
+ RTYPE ("init-wl-delta",expand->init_wl_delta,1.0);
+ EETYPE("wl-oneovert",expand->bWLoneovert,yesno_names);
+
+ *ninp_p = ninp;
+ *inp_p = inp;
+
+ return;
+}
+
void get_ir(const char *mdparin,const char *mdparout,
t_inputrec *ir,t_gromppopts *opts,
warninp_t wi)
const char *tmp;
int i,j,m,ninp;
char warn_buf[STRLEN];
-
+ t_lambda *fep = ir->fepvals;
+ t_expanded *expand = ir->expandedvals;
+
inp = read_inpfile(mdparin, &ninp, NULL, wi);
snew(dumstr[0],STRLEN);
snew(dumstr[1],STRLEN);
+ /* remove the following deprecated commands */
REM_TYPE("title");
REM_TYPE("cpp");
REM_TYPE("domain-decomposition");
- REPL_TYPE("unconstrained-start","continuation");
+ REM_TYPE("andersen-seed");
+ REM_TYPE("dihre");
+ REM_TYPE("dihre-fc");
REM_TYPE("dihre-tau");
REM_TYPE("nstdihreout");
REM_TYPE("nstcheckpoint");
+ /* replace the following commands with the clearer new versions*/
+ REPL_TYPE("unconstrained-start","continuation");
+ REPL_TYPE("foreign-lambda","fep-lambdas");
+
CCTYPE ("VARIOUS PREPROCESSING OPTIONS");
CTYPE ("Preprocessor information: use cpp syntax.");
CTYPE ("e.g.: -I/home/joe/doe -I/home/mary/roe");
RTYPE ("dt", ir->delta_t, 0.001);
STEPTYPE ("nsteps", ir->nsteps, 0);
CTYPE ("For exact run continuation or redoing part of a run");
- STEPTYPE ("init_step",ir->init_step, 0);
+ STEPTYPE ("init-step",ir->init_step, 0);
CTYPE ("Part index is updated automatically on checkpointing (keeps files separate)");
- ITYPE ("simulation_part", ir->simulation_part, 1);
+ ITYPE ("simulation-part", ir->simulation_part, 1);
CTYPE ("mode for center of mass motion removal");
EETYPE("comm-mode", ir->comm_mode, ecm_names);
CTYPE ("number of steps for center of mass motion removal");
CTYPE ("Force tolerance and initial step-size");
RTYPE ("emtol", ir->em_tol, 10.0);
RTYPE ("emstep", ir->em_stepsize,0.01);
- CTYPE ("Max number of iterations in relax_shells");
+ CTYPE ("Max number of iterations in relax-shells");
ITYPE ("niter", ir->niter, 20);
CTYPE ("Step size (ps^2) for minimization of flexible constraints");
RTYPE ("fcstep", ir->fc_stepsize, 0);
/* Output options */
CCTYPE ("OUTPUT CONTROL OPTIONS");
CTYPE ("Output frequency for coords (x), velocities (v) and forces (f)");
- ITYPE ("nstxout", ir->nstxout, 100);
- ITYPE ("nstvout", ir->nstvout, 100);
+ ITYPE ("nstxout", ir->nstxout, 0);
+ ITYPE ("nstvout", ir->nstvout, 0);
ITYPE ("nstfout", ir->nstfout, 0);
ir->nstcheckpoint = 1000;
CTYPE ("Output frequency for energies to log file and energy file");
- ITYPE ("nstlog", ir->nstlog, 100);
+ ITYPE ("nstlog", ir->nstlog, 1000);
ITYPE ("nstcalcenergy",ir->nstcalcenergy, -1);
ITYPE ("nstenergy", ir->nstenergy, 100);
CTYPE ("Output frequency and precision for .xtc file");
ir->ndelta = 2;
CTYPE ("Periodic boundary conditions: xyz, no, xy");
EETYPE("pbc", ir->ePBC, epbc_names);
- EETYPE("periodic_molecules", ir->bPeriodicMols, yesno_names);
+ EETYPE("periodic-molecules", ir->bPeriodicMols, yesno_names);
CTYPE ("nblist cut-off");
- RTYPE ("rlist", ir->rlist, 1.0);
+ RTYPE ("rlist", ir->rlist, -1);
CTYPE ("long-range cut-off for switched potentials");
RTYPE ("rlistlong", ir->rlistlong, -1);
EETYPE("coulombtype", ir->coulombtype, eel_names);
CTYPE ("cut-off lengths");
RTYPE ("rcoulomb-switch", ir->rcoulomb_switch, 0.0);
- RTYPE ("rcoulomb", ir->rcoulomb, 1.0);
+ RTYPE ("rcoulomb", ir->rcoulomb, -1);
CTYPE ("Relative dielectric constant for the medium and the reaction field");
- RTYPE ("epsilon_r", ir->epsilon_r, 1.0);
- RTYPE ("epsilon_rf", ir->epsilon_rf, 1.0);
+ RTYPE ("epsilon-r", ir->epsilon_r, 1.0);
+ RTYPE ("epsilon-rf", ir->epsilon_rf, 0.0);
CTYPE ("Method for doing Van der Waals");
EETYPE("vdw-type", ir->vdwtype, evdw_names);
CTYPE ("cut-off lengths");
RTYPE ("rvdw-switch", ir->rvdw_switch, 0.0);
- RTYPE ("rvdw", ir->rvdw, 1.0);
+ RTYPE ("rvdw", ir->rvdw, -1);
CTYPE ("Apply long range dispersion corrections for Energy and Pressure");
EETYPE("DispCorr", ir->eDispCorr, edispc_names);
CTYPE ("Extension of the potential lookup tables beyond the cut-off");
RTYPE ("table-extension", ir->tabext, 1.0);
CTYPE ("Seperate tables between energy group pairs");
- STYPE ("energygrp_table", egptable, NULL);
+ STYPE ("energygrp-table", egptable, NULL);
CTYPE ("Spacing for the PME/PPPM FFT grid");
RTYPE ("fourierspacing", opts->fourierspacing,0.12);
CTYPE ("FFT grid size, when a value is 0 fourierspacing will be used");
- ITYPE ("fourier_nx", ir->nkx, 0);
- ITYPE ("fourier_ny", ir->nky, 0);
- ITYPE ("fourier_nz", ir->nkz, 0);
+ ITYPE ("fourier-nx", ir->nkx, 0);
+ ITYPE ("fourier-ny", ir->nky, 0);
+ ITYPE ("fourier-nz", ir->nkz, 0);
CTYPE ("EWALD/PME/PPPM parameters");
- ITYPE ("pme_order", ir->pme_order, 4);
- RTYPE ("ewald_rtol", ir->ewald_rtol, 0.00001);
- EETYPE("ewald_geometry", ir->ewald_geometry, eewg_names);
- RTYPE ("epsilon_surface", ir->epsilon_surface, 0.0);
- EETYPE("optimize_fft",ir->bOptFFT, yesno_names);
+ ITYPE ("pme-order", ir->pme_order, 4);
+ RTYPE ("ewald-rtol", ir->ewald_rtol, 0.00001);
+ EETYPE("ewald-geometry", ir->ewald_geometry, eewg_names);
+ RTYPE ("epsilon-surface", ir->epsilon_surface, 0.0);
+ EETYPE("optimize-fft",ir->bOptFFT, yesno_names);
CCTYPE("IMPLICIT SOLVENT ALGORITHM");
- EETYPE("implicit_solvent", ir->implicit_solvent, eis_names);
+ EETYPE("implicit-solvent", ir->implicit_solvent, eis_names);
CCTYPE ("GENERALIZED BORN ELECTROSTATICS");
CTYPE ("Algorithm for calculating Born radii");
- EETYPE("gb_algorithm", ir->gb_algorithm, egb_names);
+ EETYPE("gb-algorithm", ir->gb_algorithm, egb_names);
CTYPE ("Frequency of calculating the Born radii inside rlist");
ITYPE ("nstgbradii", ir->nstgbradii, 1);
CTYPE ("Cutoff for Born radii calculation; the contribution from atoms");
CTYPE ("between rlist and rgbradii is updated every nstlist steps");
RTYPE ("rgbradii", ir->rgbradii, 1.0);
CTYPE ("Dielectric coefficient of the implicit solvent");
- RTYPE ("gb_epsilon_solvent",ir->gb_epsilon_solvent, 80.0);
+ RTYPE ("gb-epsilon-solvent",ir->gb_epsilon_solvent, 80.0);
CTYPE ("Salt concentration in M for Generalized Born models");
- RTYPE ("gb_saltconc", ir->gb_saltconc, 0.0);
+ RTYPE ("gb-saltconc", ir->gb_saltconc, 0.0);
CTYPE ("Scaling factors used in the OBC GB model. Default values are OBC(II)");
- RTYPE ("gb_obc_alpha", ir->gb_obc_alpha, 1.0);
- RTYPE ("gb_obc_beta", ir->gb_obc_beta, 0.8);
- RTYPE ("gb_obc_gamma", ir->gb_obc_gamma, 4.85);
- RTYPE ("gb_dielectric_offset", ir->gb_dielectric_offset, 0.009);
- EETYPE("sa_algorithm", ir->sa_algorithm, esa_names);
+ RTYPE ("gb-obc-alpha", ir->gb_obc_alpha, 1.0);
+ RTYPE ("gb-obc-beta", ir->gb_obc_beta, 0.8);
+ RTYPE ("gb-obc-gamma", ir->gb_obc_gamma, 4.85);
+ RTYPE ("gb-dielectric-offset", ir->gb_dielectric_offset, 0.009);
+ EETYPE("sa-algorithm", ir->sa_algorithm, esa_names);
CTYPE ("Surface tension (kJ/mol/nm^2) for the SA (nonpolar surface) part of GBSA");
CTYPE ("The value -1 will set default value for Still/HCT/OBC GB-models.");
- RTYPE ("sa_surface_tension", ir->sa_surface_tension, -1);
+ RTYPE ("sa-surface-tension", ir->sa_surface_tension, -1);
/* Coupling stuff */
CCTYPE ("OPTIONS FOR WEAK COUPLING ALGORITHMS");
EETYPE("tcoupl", ir->etc, etcoupl_names);
ITYPE ("nsttcouple", ir->nsttcouple, -1);
ITYPE("nh-chain-length", ir->opts.nhchainlength, NHCHAINLENGTH);
+ EETYPE("print-nose-hoover-chain-variables", ir->bPrintNHChains, yesno_names);
CTYPE ("Groups to couple separately");
STYPE ("tc-grps", tcgrps, NULL);
CTYPE ("Time constant (ps) and reference temperature (K)");
STYPE ("tau-t", tau_t, NULL);
STYPE ("ref-t", ref_t, NULL);
- CTYPE ("Pressure coupling");
- EETYPE("Pcoupl", ir->epc, epcoupl_names);
- EETYPE("Pcoupltype", ir->epct, epcoupltype_names);
+ CTYPE ("pressure coupling");
+ EETYPE("pcoupl", ir->epc, epcoupl_names);
+ EETYPE("pcoupltype", ir->epct, epcoupltype_names);
ITYPE ("nstpcouple", ir->nstpcouple, -1);
CTYPE ("Time constant (ps), compressibility (1/bar) and reference P (bar)");
RTYPE ("tau-p", ir->tau_p, 1.0);
STYPE ("compressibility", dumstr[0], NULL);
STYPE ("ref-p", dumstr[1], NULL);
CTYPE ("Scaling of reference coordinates, No, All or COM");
- EETYPE ("refcoord_scaling",ir->refcoord_scaling,erefscaling_names);
-
- CTYPE ("Random seed for Andersen thermostat");
- ITYPE ("andersen_seed", ir->andersen_seed, 815131);
+ EETYPE ("refcoord-scaling",ir->refcoord_scaling,erefscaling_names);
/* QMMM */
CCTYPE ("OPTIONS FOR QMMM calculations");
CTYPE ("Type of annealing for each temperature group (no/single/periodic)");
STYPE ("annealing", anneal, NULL);
CTYPE ("Number of time points to use for specifying annealing in each group");
- STYPE ("annealing_npoints", anneal_npoints, NULL);
+ STYPE ("annealing-npoints", anneal_npoints, NULL);
CTYPE ("List of times at the annealing points for each group");
- STYPE ("annealing_time", anneal_time, NULL);
+ STYPE ("annealing-time", anneal_time, NULL);
CTYPE ("Temp. at each annealing point, for each group.");
- STYPE ("annealing_temp", anneal_temp, NULL);
+ STYPE ("annealing-temp", anneal_temp, NULL);
/* Startup run */
CCTYPE ("GENERATE VELOCITIES FOR STARTUP RUN");
/* Energy group exclusions */
CCTYPE ("ENERGY GROUP EXCLUSIONS");
CTYPE ("Pairs of energy groups for which all non-bonded interactions are excluded");
- STYPE ("energygrp_excl", egpexcl, NULL);
+ STYPE ("energygrp-excl", egpexcl, NULL);
/* Walls */
CCTYPE ("WALLS");
CTYPE ("Number of walls, type, atom types, densities and box-z scale factor for Ewald");
ITYPE ("nwall", ir->nwall, 0);
- EETYPE("wall_type", ir->wall_type, ewt_names);
- RTYPE ("wall_r_linpot", ir->wall_r_linpot, -1);
- STYPE ("wall_atomtype", wall_atomtype, NULL);
- STYPE ("wall_density", wall_density, NULL);
- RTYPE ("wall_ewald_zfac", ir->wall_ewald_zfac, 3);
+ EETYPE("wall-type", ir->wall_type, ewt_names);
+ RTYPE ("wall-r-linpot", ir->wall_r_linpot, -1);
+ STYPE ("wall-atomtype", wall_atomtype, NULL);
+ STYPE ("wall-density", wall_density, NULL);
+ RTYPE ("wall-ewald-zfac", ir->wall_ewald_zfac, 3);
/* COM pulling */
CCTYPE("COM PULLING");
- CTYPE("Pull type: no, umbrella, constraint or constant_force");
+ CTYPE("Pull type: no, umbrella, constraint or constant-force");
EETYPE("pull", ir->ePull, epull_names);
if (ir->ePull != epullNO) {
snew(ir->pull,1);
pull_grp = read_pullparams(&ninp,&inp,ir->pull,&opts->pull_start,wi);
}
+
+ /* Enforced rotation */
+ CCTYPE("ENFORCED ROTATION");
+ CTYPE("Enforced rotation: No or Yes");
+ EETYPE("rotation", ir->bRot, yesno_names);
+ if (ir->bRot) {
+ snew(ir->rot,1);
+ rot_grp = read_rotparams(&ninp,&inp,ir->rot,wi);
+ }
/* Refinement */
CCTYPE("NMR refinement stuff");
STYPE ("orire-fitgrp",orirefitgrp, NULL);
CTYPE ("Output frequency for trace(SD) and S to energy file");
ITYPE ("nstorireout", ir->nstorireout, 100);
- CTYPE ("Dihedral angle restraints: No or Yes");
- EETYPE("dihre", opts->bDihre, yesno_names);
- RTYPE ("dihre-fc", ir->dihre_fc, 1000.0);
-
- /* Free energy stuff */
- CCTYPE ("Free energy control stuff");
- EETYPE("free-energy", ir->efep, efep_names);
- RTYPE ("init-lambda", ir->init_lambda,0.0);
- RTYPE ("delta-lambda",ir->delta_lambda,0.0);
- STYPE ("foreign_lambda", foreign_lambda, NULL);
- RTYPE ("sc-alpha",ir->sc_alpha,0.0);
- ITYPE ("sc-power",ir->sc_power,0);
- RTYPE ("sc-sigma",ir->sc_sigma,0.3);
- ITYPE ("nstdhdl", ir->nstdhdl, 10);
- EETYPE("separate-dhdl-file", ir->separate_dhdl_file,
- separate_dhdl_file_names);
- EETYPE("dhdl-derivatives", ir->dhdl_derivatives, dhdl_derivatives_names);
- ITYPE ("dh_hist_size", ir->dh_hist_size, 0);
- RTYPE ("dh_hist_spacing", ir->dh_hist_spacing, 0.1);
+
+ /* free energy variables */
+ CCTYPE ("Free energy variables");
+ EETYPE("free-energy", ir->efep, efep_names);
STYPE ("couple-moltype", couple_moltype, NULL);
EETYPE("couple-lambda0", opts->couple_lam0, couple_lam);
EETYPE("couple-lambda1", opts->couple_lam1, couple_lam);
EETYPE("couple-intramol", opts->bCoupleIntra, yesno_names);
+ RTYPE ("init-lambda", fep->init_lambda,-1); /* start with -1 so
+ we can recognize if
+ it was not entered */
+ ITYPE ("init-lambda-state", fep->init_fep_state,0);
+ RTYPE ("delta-lambda",fep->delta_lambda,0.0);
+ ITYPE ("nstdhdl",fep->nstdhdl, 10);
+ STYPE ("fep-lambdas", fep_lambda[efptFEP], NULL);
+ STYPE ("mass-lambdas", fep_lambda[efptMASS], NULL);
+ STYPE ("coul-lambdas", fep_lambda[efptCOUL], NULL);
+ STYPE ("vdw-lambdas", fep_lambda[efptVDW], NULL);
+ STYPE ("bonded-lambdas", fep_lambda[efptBONDED], NULL);
+ STYPE ("restraint-lambdas", fep_lambda[efptRESTRAINT], NULL);
+ STYPE ("temperature-lambdas", fep_lambda[efptTEMPERATURE], NULL);
+ STYPE ("init-lambda-weights",lambda_weights,NULL);
+ EETYPE("dhdl-print-energy", fep->bPrintEnergy, yesno_names);
+ RTYPE ("sc-alpha",fep->sc_alpha,0.0);
+ ITYPE ("sc-power",fep->sc_power,1);
+ RTYPE ("sc-r-power",fep->sc_r_power,6.0);
+ RTYPE ("sc-sigma",fep->sc_sigma,0.3);
+ EETYPE("sc-coul",fep->bScCoul,yesno_names);
+ ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
+ RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
+ EETYPE("separate-dhdl-file", fep->separate_dhdl_file,
+ separate_dhdl_file_names);
+ EETYPE("dhdl-derivatives", fep->dhdl_derivatives, dhdl_derivatives_names);
+ ITYPE ("dh_hist_size", fep->dh_hist_size, 0);
+ RTYPE ("dh_hist_spacing", fep->dh_hist_spacing, 0.1);
+
/* Non-equilibrium MD stuff */
CCTYPE("Non-equilibrium MD stuff");
STYPE ("acc-grps", accgrps, NULL);
RTYPE ("cos-acceleration", ir->cos_accel, 0);
STYPE ("deform", deform, NULL);
+ /* simulated tempering variables */
+ CCTYPE("simulated tempering variables");
+ EETYPE("simulated-tempering",ir->bSimTemp,yesno_names);
+ EETYPE("simulated-tempering-scaling",ir->simtempvals->eSimTempScale,esimtemp_names);
+ RTYPE("sim-temp-low",ir->simtempvals->simtemp_low,300.0);
+ RTYPE("sim-temp-high",ir->simtempvals->simtemp_high,300.0);
+
+ /* expanded ensemble variables */
+ if (ir->efep==efepEXPANDED || ir->bSimTemp)
+ {
+ read_expandedparams(&ninp,&inp,expand,wi);
+ }
+
/* Electric fields */
CCTYPE("Electric fields");
CTYPE ("Format is number of terms (int) and for all terms an amplitude (real)");
STYPE ("E-z", efield_z, NULL);
STYPE ("E-zt", efield_zt, NULL);
+ /* AdResS defined thingies */
+ CCTYPE ("AdResS parameters");
+ EETYPE("adress", ir->bAdress, yesno_names);
+ if (ir->bAdress) {
+ snew(ir->adress,1);
+ read_adressparams(&ninp,&inp,ir->adress,wi);
+ }
+
/* User defined thingies */
CCTYPE ("User defined thingies");
STYPE ("user1-grps", user1, NULL);
ir->nstcomm = 0;
opts->couple_moltype = NULL;
- if (strlen(couple_moltype) > 0) {
- if (ir->efep != efepNO) {
- opts->couple_moltype = strdup(couple_moltype);
- if (opts->couple_lam0 == opts->couple_lam1)
- warning(wi,"The lambda=0 and lambda=1 states for coupling are identical");
- if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
- opts->couple_lam1 == ecouplamNONE)) {
- warning(wi,"For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
+ if (strlen(couple_moltype) > 0)
+ {
+ if (ir->efep != efepNO)
+ {
+ opts->couple_moltype = strdup(couple_moltype);
+ if (opts->couple_lam0 == opts->couple_lam1)
+ {
+ warning(wi,"The lambda=0 and lambda=1 states for coupling are identical");
+ }
+ if (ir->eI == eiMD && (opts->couple_lam0 == ecouplamNONE ||
+ opts->couple_lam1 == ecouplamNONE))
+ {
+ warning(wi,"For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used");
+ }
+ }
+ else
+ {
+ warning(wi,"Can not couple a molecule with free_energy = no");
}
- } else {
- warning(wi,"Can not couple a molecule with free_energy = no");
- }
+ }
+ /* FREE ENERGY AND EXPANDED ENSEMBLE OPTIONS */
+ if (ir->efep != efepNO) {
+ if (fep->delta_lambda > 0) {
+ ir->efep = efepSLOWGROWTH;
+ }
+ }
+
+ if (ir->bSimTemp) {
+ fep->bPrintEnergy = TRUE;
+ /* always print out the energy to dhdl if we are doing expanded ensemble, since we need the total energy
+ if the temperature is changing. */
}
+ if ((ir->efep != efepNO) || ir->bSimTemp)
+ {
+ ir->bExpanded = FALSE;
+ if ((ir->efep == efepEXPANDED) || ir->bSimTemp)
+ {
+ ir->bExpanded = TRUE;
+ }
+ do_fep_params(ir,fep_lambda,lambda_weights);
+ if (ir->bSimTemp) { /* done after fep params */
+ do_simtemp_params(ir);
+ }
+ }
+ else
+ {
+ ir->fepvals->n_lambda = 0;
+ }
+
+ /* WALL PARAMETERS */
+
do_wall_params(ir,wall_atomtype,wall_density,opts);
+
+ /* ORIENTATION RESTRAINT PARAMETERS */
if (opts->bOrire && str_nelem(orirefitgrp,MAXPTR,NULL)!=1) {
warning_error(wi,"ERROR: Need one orientation restraint fit group\n");
}
+ /* DEFORMATION PARAMETERS */
+
clear_mat(ir->deform);
for(i=0; i<6; i++)
- dumdub[0][i] = 0;
+ {
+ dumdub[0][i] = 0;
+ }
m = sscanf(deform,"%lf %lf %lf %lf %lf %lf",
&(dumdub[0][0]),&(dumdub[0][1]),&(dumdub[0][2]),
&(dumdub[0][3]),&(dumdub[0][4]),&(dumdub[0][5]));
for(i=0; i<3; i++)
- ir->deform[i][i] = dumdub[0][i];
+ {
+ ir->deform[i][i] = dumdub[0][i];
+ }
ir->deform[YY][XX] = dumdub[0][3];
ir->deform[ZZ][XX] = dumdub[0][4];
ir->deform[ZZ][YY] = dumdub[0][5];
}
}
- if (ir->efep != efepNO) {
- parse_n_double(foreign_lambda,&ir->n_flambda,&ir->flambda);
- if (ir->n_flambda > 0 && ir->rlist < max(ir->rvdw,ir->rcoulomb)) {
- warning_note(wi,"For foreign lambda free energy differences it is assumed that the soft-core interactions have no effect beyond the neighborlist cut-off");
- }
- } else {
- ir->n_flambda = 0;
- }
-
sfree(dumstr[0]);
sfree(dumstr[1]);
}
ia = molt->ilist[F_SETTLE].iatoms;
for(i=0; i<molt->ilist[F_SETTLE].nr; ) {
/* Subtract 1 dof from every atom in the SETTLE */
- for(ai=as+ia[1]; ai<as+ia[1]+3; ai++) {
+ for(j=0; j<3; j++) {
+ ai = as + ia[1+j];
imin = min(2,nrdf2[ai]);
nrdf2[ai] -= imin;
nrdf_tc [ggrpnr(groups,egcTC ,ai)] -= 0.5*imin;
nrdf_vcm[ggrpnr(groups,egcVCM,ai)] -= 0.5*imin;
}
- ia += 2;
- i += 2;
+ ia += 4;
+ i += 4;
}
as += molt->atoms.nr;
}
nref_t = str_nelem(ref_t,MAXPTR,ptr2);
ntcg = str_nelem(tcgrps,MAXPTR,ptr3);
if ((ntau_t != ntcg) || (nref_t != ntcg)) {
- gmx_fatal(FARGS,"Invalid T coupling input: %d groups, %d ref_t values and "
- "%d tau_t values",ntcg,nref_t,ntau_t);
+ gmx_fatal(FARGS,"Invalid T coupling input: %d groups, %d ref-t values and "
+ "%d tau-t values",ntcg,nref_t,ntau_t);
}
bSetTCpar = (ir->etc || EI_SD(ir->eI) || ir->eI==eiBD || EI_TPI(ir->eI));
snew(ir->opts.tau_t,nr);
snew(ir->opts.ref_t,nr);
if (ir->eI==eiBD && ir->bd_fric==0) {
- fprintf(stderr,"bd_fric=0, so tau_t will be used as the inverse friction constant(s)\n");
+ fprintf(stderr,"bd-fric=0, so tau-t will be used as the inverse friction constant(s)\n");
}
if (bSetTCpar)
{
if (nr != nref_t)
{
- gmx_fatal(FARGS,"Not enough ref_t and tau_t values!");
+ gmx_fatal(FARGS,"Not enough ref-t and tau-t values!");
}
tau_min = 1e20;
ir->opts.tau_t[i] = strtod(ptr1[i],NULL);
if ((ir->eI == eiBD || ir->eI == eiSD2) && ir->opts.tau_t[i] <= 0)
{
- sprintf(warn_buf,"With integrator %s tau_t should be larger than 0",ei_names[ir->eI]);
+ sprintf(warn_buf,"With integrator %s tau-t should be larger than 0",ei_names[ir->eI]);
warning_error(wi,warn_buf);
}
if ((ir->etc == etcVRESCALE && ir->opts.tau_t[i] >= 0) ||
{
ir->nsttcouple = ir_optimal_nsttcouple(ir);
}
+
if (EI_VV(ir->eI))
{
+ if ((ir->etc==etcNOSEHOOVER) && (ir->epc==epcBERENDSEN)) {
+ gmx_fatal(FARGS,"Cannot do Nose-Hoover temperature with Berendsen pressure control with md-vv; use either vrescale temperature with berendsen pressure or Nose-Hoover temperature with MTTK pressure");
+ }
if ((ir->epc==epcMTTK) && (ir->etc>etcNO))
{
int mincouple;
}
ir->nstpcouple = mincouple;
ir->nsttcouple = mincouple;
- warning_note(wi,"for current Trotter decomposition methods with vv, nsttcouple and nstpcouple must be equal. Both have been reset to min(nsttcouple,nstpcouple)");
+ sprintf(warn_buf,"for current Trotter decomposition methods with vv, nsttcouple and nstpcouple must be equal. Both have been reset to min(nsttcouple,nstpcouple) = %d",mincouple);
+ warning_note(wi,warn_buf);
+ }
+ }
+ /* velocity verlet with averaged kinetic energy KE = 0.5*(v(t+1/2) - v(t-1/2)) is implemented
+ primarily for testing purposes, and does not work with temperature coupling other than 1 */
+
+ if (ETC_ANDERSEN(ir->etc)) {
+ if (ir->nsttcouple != 1) {
+ ir->nsttcouple = 1;
+ sprintf(warn_buf,"Andersen temperature control methods assume nsttcouple = 1; there is no need for larger nsttcouple > 1, since no global parameters are computed. nsttcouple has been reset to 1");
+ warning_note(wi,warn_buf);
}
}
nstcmin = tcouple_min_integration_steps(ir->etc);
{
if (tau_min/(ir->delta_t*ir->nsttcouple) < nstcmin)
{
- sprintf(warn_buf,"For proper integration of the %s thermostat, tau_t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
+ sprintf(warn_buf,"For proper integration of the %s thermostat, tau-t (%g) should be at least %d times larger than nsttcouple*dt (%g)",
ETCOUPLTYPE(ir->etc),
tau_min,nstcmin,
ir->nsttcouple*ir->delta_t);
ir->opts.ref_t[i] = strtod(ptr2[i],NULL);
if (ir->opts.ref_t[i] < 0)
{
- gmx_fatal(FARGS,"ref_t for group %d negative",i);
+ gmx_fatal(FARGS,"ref-t for group %d negative",i);
}
}
+ /* set the lambda mc temperature to the md integrator temperature (which should be defined
+ if we are in this conditional) if mc_temp is negative */
+ if (ir->expandedvals->mc_temp < 0)
+ {
+ ir->expandedvals->mc_temp = ir->opts.ref_t[0]; /*for now, set to the first reft */
+ }
}
-
+
/* Simulated annealing for each group. There are nr groups */
nSA = str_nelem(anneal,MAXPTR,ptr1);
if (nSA == 1 && (ptr1[0][0]=='n' || ptr1[0][0]=='N'))
/* Read the other fields too */
nSA_points = str_nelem(anneal_npoints,MAXPTR,ptr1);
if(nSA_points!=nSA)
- gmx_fatal(FARGS,"Found %d annealing_npoints values for %d groups\n",nSA_points,nSA);
+ gmx_fatal(FARGS,"Found %d annealing-npoints values for %d groups\n",nSA_points,nSA);
for(k=0,i=0;i<nr;i++) {
ir->opts.anneal_npoints[i]=strtol(ptr1[i],NULL,10);
if(ir->opts.anneal_npoints[i]==1)
nSA_time = str_nelem(anneal_time,MAXPTR,ptr1);
if(nSA_time!=k)
- gmx_fatal(FARGS,"Found %d annealing_time values, wanter %d\n",nSA_time,k);
+ gmx_fatal(FARGS,"Found %d annealing-time values, wanter %d\n",nSA_time,k);
nSA_temp = str_nelem(anneal_temp,MAXPTR,ptr2);
if(nSA_temp!=k)
- gmx_fatal(FARGS,"Found %d annealing_temp values, wanted %d\n",nSA_temp,k);
+ gmx_fatal(FARGS,"Found %d annealing-temp values, wanted %d\n",nSA_temp,k);
for(i=0,k=0;i<nr;i++) {
if (ir->ePull != epullNO) {
make_pull_groups(ir->pull,pull_grp,grps,gnames);
}
+
+ if (ir->bRot) {
+ make_rotation_groups(ir->rot,rot_grp,grps,gnames);
+ }
nacc = str_nelem(acc,MAXPTR,ptr1);
nacg = str_nelem(accgrps,MAXPTR,ptr2);
nr = groups->grps[egcENER].nr;
snew(ir->opts.egp_flags,nr*nr);
- bExcl = do_egp_flag(ir,groups,"energygrp_excl",egpexcl,EGP_EXCL);
+ bExcl = do_egp_flag(ir,groups,"energygrp-excl",egpexcl,EGP_EXCL);
if (bExcl && EEL_FULL(ir->coulombtype))
warning(wi,"Can not exclude the lattice Coulomb energy between energy groups");
- bTable = do_egp_flag(ir,groups,"energygrp_table",egptable,EGP_TABLE);
+ bTable = do_egp_flag(ir,groups,"energygrp-table",egptable,EGP_TABLE);
if (bTable && !(ir->vdwtype == evdwUSER) &&
!(ir->coulombtype == eelUSER) && !(ir->coulombtype == eelPMEUSER) &&
!(ir->coulombtype == eelPMEUSERSWITCH))
decode_cos(efield_yt,&(ir->et[YY]),TRUE);
decode_cos(efield_z,&(ir->ex[ZZ]),FALSE);
decode_cos(efield_zt,&(ir->et[ZZ]),TRUE);
-
+
+ if (ir->bAdress)
+ do_adress_index(ir->adress,groups,gnames,&(ir->opts),wi);
+
for(i=0; (i<grps->nr); i++)
sfree(gnames[i]);
sfree(gnames);
}
else {
sprintf(err_buf,"When using coulombtype = %s"
- " ref_t for temperature coupling should be > 0",
+ " ref-t for temperature coupling should be > 0",
eel_names[eelGRF]);
CHECK((ir->coulombtype == eelGRF) && (ir->opts.ref_t[0] <= 0));
}
-
- if (ir->eI == eiSD1) {
- gdt_max = 0;
- for(i=0; (i<ir->opts.ngtc); i++)
- gdt_max = max(gdt_max,ir->delta_t/ir->opts.tau_t[i]);
- if (0.5*gdt_max > 0.0015) {
- sprintf(warn_buf,"The relative error with integrator %s is 0.5*delta_t/tau_t = %g, you might want to switch to integrator %s\n",
- ei_names[ir->eI],0.5*gdt_max,ei_names[eiSD2]);
- warning_note(wi,warn_buf);
- }
- }
+ if (ir->eI == eiSD1 &&
+ (gmx_mtop_ftype_count(sys,F_CONSTR) > 0 ||
+ gmx_mtop_ftype_count(sys,F_SETTLE) > 0))
+ {
+ sprintf(warn_buf,"With constraints integrator %s is less accurate, consider using %s instead",ei_names[ir->eI],ei_names[eiSD2]);
+ warning_note(wi,warn_buf);
+ }
+
bAcc = FALSE;
for(i=0; (i<sys->groups.grps[egcACC].nr); i++) {
for(m=0; (m<DIM); m++) {
sfree(mgrp);
}
- if (ir->efep != efepNO && ir->sc_alpha != 0 &&
+ if (ir->efep != efepNO && ir->fepvals->sc_alpha != 0 &&
!gmx_within_tol(sys->ffparams.reppow,12.0,10*GMX_DOUBLE_EPS)) {
gmx_fatal(FARGS,"Soft-core interactions are only supported with VdW repulsion power 12");
}
if (bConstr && ir->eConstrAlg == econtSHAKE) {
if (ir->shake_tol <= 0.0) {
- sprintf(warn_buf,"ERROR: shake_tol must be > 0 instead of %g\n",
+ sprintf(warn_buf,"ERROR: shake-tol must be > 0 instead of %g\n",
ir->shake_tol);
warning_error(wi,warn_buf);
}
}
if ((ir->eI == eiCG || ir->eI == eiLBFGS) && (ir->nProjOrder<8)) {
- sprintf(warn_buf,"For accurate %s with LINCS constraints, lincs_order should be 8 or more.",ei_names[ir->eI]);
+ sprintf(warn_buf,"For accurate %s with LINCS constraints, lincs-order should be 8 or more.",ei_names[ir->eI]);
warning_note(wi,warn_buf);
}
if (ir->epc==epcMTTK) {
#include "domdec.h"
#include "partdec.h"
+#define PROBABILITYCUTOFF 100
+/* we don't bother evaluating if events are more rare than exp(-100) = 3.7x10^-44 */
+
+enum { ereTEMP, ereLAMBDA, ereENDSINGLE ,ereTL, ereNR };
+const char *erename[ereNR] = { "temperature", "lambda", "end_single_marker", "temperature and lambda"};
+
typedef struct gmx_repl_ex
{
int repl;
int nrepl;
real temp;
int type;
- real *q;
+ real **q;
gmx_bool bNPT;
real *pres;
int *ind;
+ int *allswaps;
int nst;
+ int nex;
int seed;
int nattempt[2];
real *prob_sum;
+ int **nmoves;
int *nexchange;
} t_gmx_repl_ex;
-enum { ereTEMP, ereLAMBDA, ereNR };
-const char *erename[ereNR] = { "temperature", "lambda" };
-
-static void repl_quantity(FILE *fplog,const gmx_multisim_t *ms,
- struct gmx_repl_ex *re,int ere,real q)
+static gmx_bool repl_quantity(FILE *fplog,const gmx_multisim_t *ms,
+ struct gmx_repl_ex *re,int ere,real q)
{
real *qall;
gmx_bool bDiff;
- int s;
+ int i,s;
snew(qall,ms->nsim);
qall[re->repl] = q;
gmx_sum_sim(ms->nsim,qall,ms);
bDiff = FALSE;
- for(s=1; s<ms->nsim; s++)
+ for (s=1; s<ms->nsim; s++)
{
if (qall[s] != qall[0])
{
- bDiff = TRUE;
+ bDiff = TRUE;
}
}
+
if (bDiff)
{
- if (re->type >= 0 && re->type < ereNR)
- {
- gmx_fatal(FARGS,"For replica exchange both %s and %s differ",
- erename[re->type],erename[ere]);
- }
/* Set the replica exchange type and quantities */
re->type = ere;
- snew(re->q,re->nrepl);
+
+ snew(re->q[ere],re->nrepl);
for(s=0; s<ms->nsim; s++)
{
- re->q[s] = qall[s];
+ re->q[ere][s] = qall[s];
}
}
-
sfree(qall);
+ return bDiff;
}
gmx_repl_ex_t init_replica_exchange(FILE *fplog,
const gmx_multisim_t *ms,
const t_state *state,
const t_inputrec *ir,
- int nst,int init_seed)
+ int nst, int nex, int init_seed)
{
real temp,pres;
int i,j,k;
struct gmx_repl_ex *re;
+ gmx_bool bTemp;
+ gmx_bool bLambda=FALSE;
fprintf(fplog,"\nInitializing Replica Exchange\n");
re->repl = ms->sim;
re->nrepl = ms->nsim;
+ snew(re->q,ereENDSINGLE);
fprintf(fplog,"Repl There are %d replicas:\n",re->nrepl);
"the number of temperature coupling groups");
check_multi_int(fplog,ms,ir->epc,"the pressure coupling");
check_multi_int(fplog,ms,ir->efep,"free energy");
+ check_multi_int(fplog,ms,ir->fepvals->n_lambda,"number of lambda states");
re->temp = ir->opts.ref_t[0];
for(i=1; (i<ir->opts.ngtc); i++)
}
re->type = -1;
- for(i=0; i<ereNR; i++)
+ bTemp = repl_quantity(fplog,ms,re,ereTEMP,re->temp);
+ if (ir->efep != efepNO)
{
- switch (i)
- {
- case ereTEMP:
- repl_quantity(fplog,ms,re,i,re->temp);
- break;
- case ereLAMBDA:
- if (ir->efep != efepNO)
- {
- repl_quantity(fplog,ms,re,i,ir->init_lambda);
- }
- break;
- default:
- gmx_incons("Unknown replica exchange quantity");
- }
+ bLambda = repl_quantity(fplog,ms,re,ereLAMBDA,(real)ir->fepvals->init_fep_state);
}
- if (re->type == -1)
+ if (re->type == -1) /* nothing was assigned */
{
gmx_fatal(FARGS,"The properties of the %d systems are all the same, there is nothing to exchange",re->nrepl);
}
+ if (bLambda && bTemp) {
+ re->type = ereTL;
+ }
- switch (re->type)
+ if (bTemp)
{
- please_cite(fplog,"Hukushima96a");
- case ereTEMP:
+ please_cite(fplog,"Sugita1999a");
if (ir->epc != epcNO)
{
re->bNPT = TRUE;
gmx_fatal(FARGS,"REMD with the %s thermostat does not produce correct potential energy distributions, consider using the %s thermostat instead",
ETCOUPLTYPE(ir->etc),ETCOUPLTYPE(etcVRESCALE));
}
- break;
- case ereLAMBDA:
- if (ir->delta_lambda != 0)
+ }
+ if (bLambda) {
+ if (ir->fepvals->delta_lambda != 0) /* check this? */
{
gmx_fatal(FARGS,"delta_lambda is not zero");
}
- break;
}
-
if (re->bNPT)
{
snew(re->pres,re->nrepl);
gmx_sum_sim(re->nrepl,re->pres,ms);
}
+ /* Make an index for increasing replica order */
+ /* only makes sense if one or the other is varying, not both!
+ if both are varying, we trust the order the person gave. */
snew(re->ind,re->nrepl);
- /* Make an index for increasing temperature order */
for(i=0; i<re->nrepl; i++)
{
re->ind[i] = i;
}
- for(i=0; i<re->nrepl; i++)
- {
- for(j=i+1; j<re->nrepl; j++)
+
+ if (re->type<ereENDSINGLE) {
+
+ for(i=0; i<re->nrepl; i++)
{
- if (re->q[re->ind[j]] < re->q[re->ind[i]])
+ for(j=i+1; j<re->nrepl; j++)
{
- k = re->ind[i];
- re->ind[i] = re->ind[j];
- re->ind[j] = k;
- }
- else if (re->q[re->ind[j]] == re->q[re->ind[i]])
- {
- gmx_fatal(FARGS,"Two replicas have identical %ss",erename[re->type]);
+ if (re->q[re->type][re->ind[j]] < re->q[re->type][re->ind[i]])
+ {
+ k = re->ind[i];
+ re->ind[i] = re->ind[j];
+ re->ind[j] = k;
+ }
+ else if (re->q[re->type][re->ind[j]] == re->q[re->type][re->ind[i]])
+ {
+ gmx_fatal(FARGS,"Two replicas have identical %ss",erename[re->type]);
+ }
}
}
}
- fprintf(fplog,"Repl ");
+
+ /* keep track of all the swaps, starting with the initial placement. */
+ snew(re->allswaps,re->nrepl);
for(i=0; i<re->nrepl; i++)
{
- fprintf(fplog," %3d ",re->ind[i]);
+ re->allswaps[i] = re->ind[i];
}
+
switch (re->type)
{
case ereTEMP:
- fprintf(fplog,"\nRepl T");
+ fprintf(fplog,"\nReplica exchange in temperature\n");
for(i=0; i<re->nrepl; i++)
{
- fprintf(fplog," %5.1f",re->q[re->ind[i]]);
+ fprintf(fplog," %5.1f",re->q[re->type][re->ind[i]]);
}
+ fprintf(fplog,"\n");
break;
case ereLAMBDA:
- fprintf(fplog,"\nRepl l");
+ fprintf(fplog,"\nReplica exchange in lambda\n");
+ for(i=0; i<re->nrepl; i++)
+ {
+ fprintf(fplog," %3d",(int)re->q[re->type][re->ind[i]]);
+ }
+ fprintf(fplog,"\n");
+ break;
+ case ereTL:
+ fprintf(fplog,"\nReplica exchange in temperature and lambda state\n");
+ for(i=0; i<re->nrepl; i++)
+ {
+ fprintf(fplog," %5.1f",re->q[ereTEMP][re->ind[i]]);
+ }
+ fprintf(fplog,"\n");
for(i=0; i<re->nrepl; i++)
{
- fprintf(fplog," %5.3f",re->q[re->ind[i]]);
+ fprintf(fplog," %5d",(int)re->q[ereLAMBDA][re->ind[i]]);
}
+ fprintf(fplog,"\n");
break;
default:
gmx_incons("Unknown replica exchange quantity");
{
if ((i > 0) && (re->pres[re->ind[i]] < re->pres[re->ind[i-1]]))
{
- gmx_fatal(FARGS,"The reference pressure decreases with increasing temperature");
+ fprintf(fplog,"\nWARNING: The reference pressures decrease with increasing temperatures\n\n");
+ fprintf(stderr,"\nWARNING: The reference pressures decrease with increasing temperatures\n\n");
}
}
}
- fprintf(fplog,"\nRepl ");
-
re->nst = nst;
if (init_seed == -1)
{
{
re->seed = init_seed;
}
- fprintf(fplog,"\nRepl exchange interval: %d\n",re->nst);
- fprintf(fplog,"\nRepl random seed: %d\n",re->seed);
+ fprintf(fplog,"\nReplica exchange interval: %d\n",re->nst);
+ fprintf(fplog,"\nReplica random seed: %d\n",re->seed);
re->nattempt[0] = 0;
re->nattempt[1] = 0;
+
snew(re->prob_sum,re->nrepl);
snew(re->nexchange,re->nrepl);
+ snew(re->nmoves,re->nrepl);
+ for (i=0;i<re->nrepl;i++)
+ {
+ snew(re->nmoves[i],re->nrepl);
+ }
+ fprintf(fplog,"Replica exchange information below: x=exchange, pr=probability\n");
- fprintf(fplog,"Repl below: x=exchange, pr=probability\n");
-
+ re->nex = nex;
return re;
}
}
}
+
+static void exchange_ints(const gmx_multisim_t *ms,int b,int *v,int n)
+{
+ int *buf;
+ int i;
+
+ if (v) {
+ snew(buf,n);
+#ifdef GMX_MPI
+ /*
+ MPI_Sendrecv(v, n*sizeof(int),MPI_BYTE,MSRANK(ms,b),0,
+ buf,n*sizeof(int),MPI_BYTE,MSRANK(ms,b),0,
+ ms->mpi_comm_masters,MPI_STATUS_IGNORE);
+ */
+ {
+ MPI_Request mpi_req;
+
+ MPI_Isend(v,n*sizeof(int),MPI_BYTE,MSRANK(ms,b),0,
+ ms->mpi_comm_masters,&mpi_req);
+ MPI_Recv(buf,n*sizeof(int),MPI_BYTE,MSRANK(ms,b),0,
+ ms->mpi_comm_masters,MPI_STATUS_IGNORE);
+ MPI_Wait(&mpi_req,MPI_STATUS_IGNORE);
+ }
+#endif
+ for(i=0; i<n; i++)
+ {
+ v[i] = buf[i];
+ }
+ sfree(buf);
+ }
+}
+
static void exchange_doubles(const gmx_multisim_t *ms,int b,double *v,int n)
{
double *buf;
}
}
+static void copy_reals(const real *s,real *d,int n)
+{
+ int i;
+
+ if (d != NULL)
+ {
+ for(i=0; i<n; i++)
+ {
+ d[i] = s[i];
+ }
+ }
+}
+
+static void copy_ints(const int *s,int *d,int n)
+{
+ int i;
+
+ if (d != NULL)
+ {
+ for(i=0; i<n; i++)
+ {
+ d[i] = s[i];
+ }
+ }
+}
+
#define scopy_rvecs(v,n) copy_rvecs(state->v,state_local->v,n);
#define scopy_doubles(v,n) copy_doubles(state->v,state_local->v,n);
+#define scopy_reals(v,n) copy_reals(state->v,state_local->v,n);
+#define scopy_ints(v,n) copy_ints(state->v,state_local->v,n);
static void copy_state_nonatomdata(t_state *state,t_state *state_local)
{
scopy_rvecs(x,state->natoms);
scopy_rvecs(v,state->natoms);
scopy_rvecs(sd_X,state->natoms);
+ copy_ints(&(state->fep_state),&(state_local->fep_state),1);
+ scopy_reals(lambda,efptNR);
}
static void scale_velocities(t_state *state,real fac)
}
}
+static void print_matrix(FILE *fplog,const char *leg,int n,int **nmoves, int *nattempt)
+{
+ int i,j,ntot;
+ float Tprint;
+
+ ntot = nattempt[0] + nattempt[1];
+
+ fprintf(fplog," Empirical Transition Matrix\n");
+ for (i=0;i<n;i++)
+ {
+ fprintf(fplog,"%8d",(i+1));
+ }
+ fprintf(fplog,"\n");
+ for (i=0;i<n;i++)
+ {
+ for (j=0;j<n;j++)
+ {
+ Tprint = 0.0;
+ if (nmoves[i][j] > 0)
+ {
+ Tprint = nmoves[i][j]/(2.0*ntot);
+ }
+ fprintf(fplog,"%8.4f",Tprint);
+ }
+ fprintf(fplog,"%3d\n",i);
+ }
+}
+
static void print_ind(FILE *fplog,const char *leg,int n,int *ind,gmx_bool *bEx)
{
int i;
fprintf(fplog,"\n");
}
+static void print_allswitchind(FILE *fplog,int n,int *ind,int *pind, int *allswaps)
+{
+ int i;
+ int *tmpswap;
+
+ snew(tmpswap,n); /* need to save the data */
+ for (i=0;i<n;i++)
+ {
+ tmpswap[i] = allswaps[i];
+ }
+ for (i=0;i<n;i++)
+ {
+ allswaps[i] = tmpswap[pind[i]];
+ }
+
+ fprintf(fplog,"\nAccepted Exchanges: ");
+ for (i=0;i<n;i++)
+ {
+ fprintf(fplog,"%d ",pind[i]);
+ }
+ fprintf(fplog,"\n");
+
+ fprintf(fplog,"Order After Exchange: ");
+ for (i=0;i<n;i++)
+ {
+ fprintf(fplog,"%d ",allswaps[i]);
+ }
+ fprintf(fplog,"\n\n");
+
+ sfree(tmpswap);
+}
+
static void print_prob(FILE *fplog,const char *leg,int n,real *prob)
{
int i;
fprintf(fplog,"\n");
}
-static int get_replica_exchange(FILE *fplog,const gmx_multisim_t *ms,
- struct gmx_repl_ex *re,real *ener,real vol,
- gmx_large_int_t step,real time)
-{
- int m,i,a,b;
- real *Epot=NULL,*Vol=NULL,*dvdl=NULL,*prob;
- real ediff=0,delta=0,dpV=0,betaA=0,betaB=0;
- gmx_bool *bEx,bPrint;
- int exchange;
+static real calc_delta(FILE *fplog, gmx_bool bPrint, struct gmx_repl_ex *re, real *Epot, real **df, real* Vol, real *beta, int a, int b, int ap, int bp) {
+
+ real ediff,dpV,delta=0;
+
+ /* Two cases; we are permuted and not. In all cases, setting ap = a and bp = b will reduce
+ to the non permuted case */
- fprintf(fplog,"Replica exchange at step " gmx_large_int_pfmt " time %g\n",step,time);
-
switch (re->type)
{
case ereTEMP:
- snew(Epot,re->nrepl);
- snew(Vol,re->nrepl);
- Epot[re->repl] = ener[F_EPOT];
- Vol[re->repl] = vol;
- gmx_sum_sim(re->nrepl,Epot,ms);
- gmx_sum_sim(re->nrepl,Vol,ms);
+ /*
+ * Okabe et. al. Chem. Phys. Lett. 335 (2001) 435-439
+ */
+ ediff = Epot[b] - Epot[a];
+ delta = -(beta[bp] - beta[ap])*ediff;
break;
case ereLAMBDA:
- snew(dvdl,re->nrepl);
- dvdl[re->repl] = ener[F_DVDL];
- gmx_sum_sim(re->nrepl,dvdl,ms);
+ /* two cases: when we are permuted, and not. */
+ /* non-permuted:
+ ediff = E_new - E_old
+ = [H_b(x_a) + H_a(x_b)] - [H_b(x_b) + H_a(x_a)]
+ = [H_b(x_a) - H_a(x_a)] + [H_a(x_b) - H_b(x_b)]
+ = df[b][a] + df[a][b] */
+ /* permuted:
+ ediff = E_new - E_old
+ = [H_bp(x_a) + H_ap(x_b)] - [H_bp(x_b) + H_ap(x_a)]
+ = [H_bp(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_bp(x_b)]
+ = [H_bp(x_a) - H_a(x_a) + H_a(x_a) - H_ap(x_a)] + [H_ap(x_b) - H_b(x_b) + H_b(x_b) - H_bp(x_b)]
+ = [H_bp(x_a) - H_a(x_a)] - [H_ap(x_a) - H_a(x_a)] + [H_ap(x_b) - H_b(x_b)] - H_bp(x_b) - H_b(x_b)]
+ = (df[bp][a] - df[ap][a]) + (df[ap][b] - df[bp][b]) */
+ ediff = (df[bp][a] - df[ap][a]) + (df[ap][b] - df[bp][b]);
+ delta = ediff*beta[a]; /* assume all same temperature in this case */
break;
+ case ereTL:
+ /* not permuted: */
+ /* delta = reduced E_new - reduced E_old
+ = [beta_b H_b(x_a) + beta_a H_a(x_b)] - [beta_b H_b(x_b) + beta_a H_a(x_a)]
+ = [beta_b H_b(x_a) - beta_a H_a(x_a)] + [beta_a H_a(x_b) - beta_b H_b(x_b)]
+ = [beta_b dH_b(x_a) + beta_b H_a(x_a) - beta_a H_a(x_a)] +
+ [beta_a dH_a(x_b) + beta_a H_b(x_b) - beta_b H_b(x_b)]
+ = [beta_b dH_b(x_a) + [beta_a dH_a(x_b) +
+ beta_b (H_a(x_a) - H_b(x_b)]) - beta_a (H_a(x_a) - H_b(x_b))
+ = beta_b dH_b(x_a) + beta_a dH_a(x_b) - (beta_b - beta_a)(H_b(x_b) - H_a(x_a) */
+ /* delta = beta[b]*df[b][a] + beta[a]*df[a][b] - (beta[b] - beta[a])*(Epot[b] - Epot[a]; */
+ /* permuted (big breath!) */
+ /* delta = reduced E_new - reduced E_old
+ = [beta_bp H_bp(x_a) + beta_ap H_ap(x_b)] - [beta_bp H_bp(x_b) + beta_ap H_ap(x_a)]
+ = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
+ = [beta_bp H_bp(x_a) - beta_ap H_ap(x_a)] + [beta_ap H_ap(x_b) - beta_bp H_bp(x_b)]
+ - beta_pb H_a(x_a) + beta_ap H_a(x_a) + beta_pb H_a(x_a) - beta_ap H_a(x_a)
+ - beta_ap H_b(x_b) + beta_bp H_b(x_b) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
+ = [(beta_bp H_bp(x_a) - beta_bp H_a(x_a)) - (beta_ap H_ap(x_a) - beta_ap H_a(x_a))] +
+ [(beta_ap H_ap(x_b) - beta_ap H_b(x_b)) - (beta_bp H_bp(x_b) - beta_bp H_b(x_b))]
+ + beta_pb H_a(x_a) - beta_ap H_a(x_a) + beta_ap H_b(x_b) - beta_bp H_b(x_b)
+ = [beta_bp (H_bp(x_a) - H_a(x_a)) - beta_ap (H_ap(x_a) - H_a(x_a))] +
+ [beta_ap (H_ap(x_b) - H_b(x_b)) - beta_bp (H_bp(x_b) - H_b(x_b))]
+ + beta_pb (H_a(x_a) - H_b(x_b)) - beta_ap (H_a(x_a) - H_b(x_b))
+ = ([beta_bp df[bp][a] - beta_ap df[ap][a]) + beta_ap df[ap][b] - beta_bp df[bp][b])
+ + (beta_pb-beta_ap)(H_a(x_a) - H_b(x_b)) */
+ delta = beta[bp]*(df[bp][a] - df[bp][b]) + beta[ap]*(df[ap][b] - df[ap][a]) - (beta[bp]-beta[ap])*(Epot[b]-Epot[a]);
+ break;
+ default:
+ gmx_incons("Unknown replica exchange quantity");
+ }
+ if (bPrint)
+ {
+ fprintf(fplog,"Repl %d <-> %d dE_term = %10.3e (kT)\n",a,b,delta);
}
- int step,real time,int *pind)
+ if (re->bNPT)
+ {
+ /* revist the calculation for 5.0. Might be some improvements. */
+ dpV = (beta[ap]*re->pres[ap]-beta[bp]*re->pres[bp])*(Vol[b]-Vol[a])/PRESFAC;
+ if (bPrint)
+ {
+ fprintf(fplog," dpV = %10.3e d = %10.3e\nb",dpV,delta + dpV);
+ }
+ delta += dpV;
+ }
+ return delta;
+}
+
+static void get_replica_exchange(FILE *fplog,const gmx_multisim_t *ms,
+ struct gmx_repl_ex *re,gmx_enerdata_t *enerd,real vol,
- fprintf(fplog,"Replica exchange at step %d time %g\n",step,time);
++ gmx_large_int_t step,real time,int *pind)
+{
+ int m,i,a,b,ap,bp,i0,i1,tmp;
+ real *Epot=NULL,*Vol=NULL,**flambda=NULL,*beta=NULL,*prob;
+ real ediff=0,delta=0,dpV=0;
+ gmx_bool *bEx,bPrint,bMultiEx;
+ gmx_bool bEpot=FALSE;
+ gmx_bool bFLambda=FALSE;
+ gmx_bool bVol=FALSE;
+
+ bMultiEx = (re->nex > 1); /* multiple exchanges at each state */
++ fprintf(fplog,"Replica exchange at step " gmx_large_int_pfmt " time %g\n",step,time);
+ snew(beta,re->nrepl);
+ if (re->bNPT)
+ {
+ bVol = TRUE;
+ snew(Vol,re->nrepl);
+ Vol[re->repl] = vol;
+ }
+
+ if ((re->type == ereTEMP || re->type == ereTL))
+ {
+ bEpot = TRUE;
+ snew(Epot,re->nrepl);
+ Epot[re->repl] = enerd->term[F_EPOT];
+ /* temperatures of different states*/
+ for (i=0;i<re->nrepl;i++)
+ {
+ beta[i] = 1.0/(re->q[ereTEMP][i]*BOLTZ);
+ }
+ }
+ else
+ {
+ for (i=0;i<re->nrepl;i++)
+ {
+ beta[i] = 1.0/(re->temp*BOLTZ); /* we have a single temperature */
+ }
+ }
+ if (re->type == ereLAMBDA || re->type == ereTL)
+ {
+ bFLambda = TRUE;
+ /* lambda differences. */
+ /* flambda[i][j] is the energy of the jth simulation in the ith Hamiltonian
+ minus the energy of the jth simulation in the jth Hamiltonian */
+ snew(flambda,re->nrepl);
+ for (i=0;i<re->nrepl;i++)
+ {
+ snew(flambda[i],re->nrepl);
+ flambda[i][re->repl] = (enerd->enerpart_lambda[(int)re->q[ereLAMBDA][i]+1]-enerd->enerpart_lambda[0]);
+ }
+ }
+
+ /* now actually do the communication */
+ if (bVol)
+ {
+ gmx_sum_sim(re->nrepl,Vol,ms);
+ }
+ if (bEpot)
+ {
+ gmx_sum_sim(re->nrepl,Epot,ms);
+ }
+ if (bFLambda)
+ {
+ for (i=0;i<re->nrepl;i++)
+ {
+ gmx_sum_sim(re->nrepl,flambda[i],ms);
+ }
+ }
snew(bEx,re->nrepl);
snew(prob,re->nrepl);
- exchange = -1;
- m = (step / re->nst) % 2;
- for(i=1; i<re->nrepl; i++)
+ /* make a duplicate set of indices for shuffling */
+ for(i=0;i<re->nrepl;i++)
+ {
+ pind[i] = re->ind[i];
+ }
+
+ if (bMultiEx)
{
- a = re->ind[i-1];
- b = re->ind[i];
- bPrint = (re->repl==a || re->repl==b);
- if (i % 2 == m)
+ /* multiple random switch exchange */
+ for (i=0;i<re->nex;i++)
{
- switch (re->type)
+ /* randomly select a pair */
+ /* find out which state it is from, and what label that state currently has */
+ i0 = (int)(re->nrepl*rando(&(re->seed)));
+ i1 = (int)(re->nrepl*rando(&(re->seed)));
+ if (i0==i1)
{
- case ereTEMP:
- /* Use equations from:
- * Okabe et. al. Chem. Phys. Lett. 335 (2001) 435-439
- */
- ediff = Epot[b] - Epot[a];
- betaA = 1.0/(re->q[a]*BOLTZ);
- betaB = 1.0/(re->q[b]*BOLTZ);
- delta = (betaA - betaB)*ediff;
- break;
- case ereLAMBDA:
- /* Here we exchange based on a linear extrapolation of dV/dlambda.
- * We would like to have the real energies
- * from foreign lambda calculations.
- */
- ediff = (dvdl[a] - dvdl[b])*(re->q[b] - re->q[a]);
- delta = ediff/(BOLTZ*re->temp);
- break;
- default:
- gmx_incons("Unknown replica exchange quantity");
- }
- if (bPrint)
- {
- fprintf(fplog,"Repl %d <-> %d dE = %10.3e",a,b,delta);
- }
- if (re->bNPT)
- {
- dpV = (betaA*re->pres[a]-betaB*re->pres[b])*(Vol[b]-Vol[a])/PRESFAC;
- if (bPrint)
- {
- fprintf(fplog," dpV = %10.3e d = %10.3e",dpV,delta + dpV);
- }
- delta += dpV;
- }
- if (bPrint)
- {
- fprintf(fplog,"\n");
+ i--;
+ continue; /* got the same pair, back up and do it again */
}
+
+ a = re->ind[i0];
+ b = re->ind[i1];
+ ap = pind[i0];
+ bp = pind[i1];
+
+ bPrint = FALSE; /* too noisy */
+ delta = calc_delta(fplog,bPrint,re,Epot,flambda,Vol,beta,a,b,ap,bp); /* calculate the energy difference */
+
+ /* we actually only use the first space, since there are actually many switches between pairs. */
+
if (delta <= 0)
{
- prob[i] = 1;
- bEx[i] = TRUE;
+ /* accepted */
+ prob[0] = 1;
+ bEx[0] = TRUE;
}
else
{
- if (delta > 100)
+ if (delta > PROBABILITYCUTOFF)
{
- prob[i] = 0;
+ prob[0] = 0;
}
else
{
- prob[i] = exp(-delta);
+ prob[0] = exp(-delta);
}
- bEx[i] = (rando(&(re->seed)) < prob[i]);
+ /* roll a number to determine if accepted */
+ bEx[0] = (rando(&(re->seed)) < prob[0]);
}
- re->prob_sum[i] += prob[i];
- if (bEx[i])
+ re->prob_sum[0] += prob[0];
+
+ if (bEx[0])
{
- if (a == re->repl)
+ /* swap the states */
+ tmp = pind[i0];
+ pind[i0] = pind[i1];
+ pind[i1] = tmp;
+ }
+ }
+ re->nattempt[0]++; /* keep track of total permutation trials here */
+ print_allswitchind(fplog,re->nrepl,re->ind,pind,re->allswaps);
+ }
+ else
+ {
+ /* standard nearest neighbor replica exchange */
+ m = (step / re->nst) % 2;
+ for(i=1; i<re->nrepl; i++)
+ {
+ a = re->ind[i-1];
+ b = re->ind[i];
+
+ bPrint = (re->repl==a || re->repl==b);
+ if (i % 2 == m)
+ {
+ delta = calc_delta(fplog,bPrint,re,Epot,flambda,Vol,beta,a,b,a,b);
+ if (delta <= 0) {
+ /* accepted */
+ prob[i] = 1;
+ bEx[i] = TRUE;
+ }
+ else
{
- exchange = b;
+ if (delta > PROBABILITYCUTOFF)
+ {
+ prob[i] = 0;
+ }
+ else
+ {
+ prob[i] = exp(-delta);
+ }
+ /* roll a number to determine if accepted */
+ bEx[i] = (rando(&(re->seed)) < prob[i]);
}
- else if (b == re->repl)
+ re->prob_sum[i] += prob[i];
+
+ if (bEx[i])
{
- exchange = a;
+ /* swap these two */
+ tmp = pind[i-1];
+ pind[i-1] = pind[i];
+ pind[i] = tmp;
}
- re->nexchange[i]++;
+ }
+ else
+ {
+ prob[i] = -1;
+ bEx[i] = FALSE;
}
}
- else
- {
- prob[i] = -1;
- bEx[i] = FALSE;
- }
+ /* print some statistics */
+ print_ind(fplog,"ex",re->nrepl,re->ind,bEx);
+ print_prob(fplog,"pr",re->nrepl,prob);
+ fprintf(fplog,"\n");
+ re->nattempt[m]++;
}
- print_ind(fplog,"ex",re->nrepl,re->ind,bEx);
- print_prob(fplog,"pr",re->nrepl,prob);
- fprintf(fplog,"\n");
+ /* record which moves were made and accepted */
+ for (i=0;i<re->nrepl;i++)
+ {
+ re->nmoves[re->ind[i]][pind[i]] +=1;
+ re->nmoves[pind[i]][re->ind[i]] +=1;
+ }
+ /* free up data */
sfree(bEx);
sfree(prob);
- sfree(Epot);
- sfree(Vol);
- sfree(dvdl);
-
- re->nattempt[m]++;
-
- return exchange;
+ sfree(beta);
+ if (re->bNPT)
+ {
+ sfree(Vol);
+ }
+ if ((re->type == ereTEMP || re->type == ereTL))
+ {
+ sfree(Epot);
+ }
+ if ((re->type == ereLAMBDA || re->type == ereTL))
+ {
+ for (i=0;i<re->nrepl;i++)
+ {
+ sfree(flambda[i]);
+ }
+ sfree(flambda);
+ }
}
static void write_debug_x(t_state *state)
}
}
+static void cyclic_decomposition(FILE *fplog, int *pind, int **cyclic, int n, int *nswap)
+{
+
+ int i,j,c,p;
+ int *incycle;
+ int maxlen = 1;
+ snew(incycle,n);
+
+ /* compute cyclic decompositions */
+ for (i=0;i<n;i++)
+ {
+ snew(cyclic[i],n);
+ for (j=0;j<n;j++)
+ {
+ cyclic[i][j] = -1;
+ }
+ }
+
+ for (i=0;i<n;i++) /* one cycle for each replica */
+ {
+ if (incycle[i])
+ {
+ cyclic[i][0] = -1;
+ continue;
+ }
+ cyclic[i][0] = i;
+ incycle[i] = TRUE;
+ c = 1;
+ p = i;
+ for (j=0;j<n;j++) /* potentially all cycles are part, but we will break first */
+ {
+ p = pind[p]; /* start permuting */
+ if (p==i)
+ {
+ cyclic[i][c] = -1;
+ if (c > maxlen)
+ {
+ maxlen = c;
+ }
+ break; /* we've reached the original element, the cycle is complete, and we marked the end. */
+ }
+ else
+ {
+ cyclic[i][c] = p; /* each permutation gives a new member of the cycle */
+ incycle[p] = TRUE;
+ c++;
+ }
+ }
+ }
+ *nswap = maxlen - 1;
+
+ if (debug)
+ {
+ for (i=0;i<n;i++)
+ {
+ fprintf(fplog,"Cycle %d:",i);
+ for (j=0;j<n;j++)
+ {
+ if (cyclic[i][j] < 0)
+ {
+ break;
+ }
+ fprintf(fplog,"%2d",cyclic[i][j]);
+ }
+ fprintf(fplog,"\n");
+ }
+ fflush(fplog);
+ }
+}
+
+static void compute_exchange_order(FILE *fplog, int **cyclic,int **order, int n, int maxswap)
+{
+ int i,j;
+
+ for (i=0;i<n;i++)
+ {
+ snew(order[i],maxswap);
+ for (j=0;j<maxswap;j++)
+ {
+ order[i][j] = -1;
+ }
+ }
+ for (j=0;j<maxswap;j++)
+ {
+ for (i=0;i<n;i++)
+ {
+ if (cyclic[i][j+1] >= 0)
+ {
+ order[cyclic[i][j+1]][j] = cyclic[i][j];
+ order[cyclic[i][j]][j] = cyclic[i][j+1];
+ }
+ }
+ for (i=0;i<n;i++)
+ {
+ if (order[i][j] < 0)
+ {
+ order[i][j] = i; /* if it's not exchanging, it should stay this round*/
+ }
+ }
+ }
+ if (debug)
+ {
+ fprintf(fplog,"Replica Exchange Order\n");
+ for (i=0;i<n;i++)
+ {
+ fprintf(fplog,"Replica %d:",i);
+ for (j=0;j<maxswap;j++)
+ {
+ if (order[i][j] < 0) break;
+ fprintf(fplog,"%2d",order[i][j]);
+ }
+ fprintf(fplog,"\n");
+ }
+ fflush(fplog);
+ }
+}
+
gmx_bool replica_exchange(FILE *fplog,const t_commrec *cr,struct gmx_repl_ex *re,
- t_state *state,real *ener,
- t_state *state_local,
- gmx_large_int_t step,real time)
+ t_state *state,gmx_enerdata_t *enerd,
- t_state *state_local,int step,real time)
++ t_state *state_local,gmx_large_int_t step,real time)
{
gmx_multisim_t *ms;
- int exchange=-1,shift;
+ int exchange=-1,shift;
+ int i,j,maxswap=0;
+ int *exchanges=NULL;
+ int **cyclic=NULL;
+ int **order=NULL;
gmx_bool bExchanged=FALSE;
-
+
ms = cr->ms;
-
if (MASTER(cr))
{
- exchange = get_replica_exchange(fplog,ms,re,ener,det(state->box),
- step,time);
- bExchanged = (exchange >= 0);
+ snew(exchanges,re->nrepl);
+ get_replica_exchange(fplog,ms,re,enerd,det(state->box),step,time,exchanges);
+ bExchanged = (exchanges[re->repl] != re->nrepl); /* only mark as exchanged if it has a shuffled index */
+ snew(cyclic,re->nrepl);
+ snew(order,re->nrepl);
+
+ /* identify the cyclic decomposition of the permutation (very easy if neighbor replica exchange) */
+ cyclic_decomposition(fplog,exchanges,cyclic,re->nrepl,&maxswap);
+
+ /* now translate the decompsition into a replica exchange order at each step */
+ compute_exchange_order(fplog,cyclic,order,re->nrepl,maxswap);
+
+ sfree(cyclic); /* don't need this anymore */
}
-
if (PAR(cr))
{
#ifdef GMX_MPI
cr->mpi_comm_mygroup);
#endif
}
-
if (bExchanged)
{
/* Exchange the states */
if (MASTER(cr))
{
- /* Exchange the global states between the master nodes */
- if (debug)
+ for (i=0;i<maxswap;i++) /* there will be only one swap cycle with standard replica exchange */
{
- fprintf(debug,"Exchanging %d with %d\n",ms->sim,exchange);
+ exchange = order[ms->sim][i];
+
+ if (exchange != ms->sim)
+ {
+ /* Exchange the global states between the master nodes */
+ if (debug)
+ {
+ fprintf(debug,"Exchanging %d with %d\n",ms->sim,exchange);
+ }
+ exchange_state(ms,exchange,state);
+ }
}
- exchange_state(ms,exchange,state);
-
- if (re->type == ereTEMP)
+ if (re->type == ereTEMP || re->type == ereTL)
{
- scale_velocities(state,sqrt(re->q[ms->sim]/re->q[exchange]));
+ scale_velocities(state,sqrt(re->q[ereTEMP][ms->sim]/re->q[ereTEMP][exchanges[ms->sim]]));
}
+ sfree(order);
}
/* With domain decomposition the global state is distributed later */
}
}
}
-
return bExchanged;
}
{
real *prob;
int i;
-
+
fprintf(fplog,"\nReplica exchange statistics\n");
- fprintf(fplog,"Repl %d attempts, %d odd, %d even\n",
- re->nattempt[0]+re->nattempt[1],re->nattempt[1],re->nattempt[0]);
- snew(prob,re->nrepl);
-
- fprintf(fplog,"Repl average probabilities:\n");
- for(i=1; i<re->nrepl; i++)
+ if (re->nex == 0)
{
- if (re->nattempt[i%2] == 0)
- {
- prob[i] = 0;
- }
- else
- {
- prob[i] = re->prob_sum[i]/re->nattempt[i%2];
- }
- }
- print_ind(fplog,"",re->nrepl,re->ind,NULL);
- print_prob(fplog,"",re->nrepl,prob);
+ fprintf(fplog,"Repl %d attempts, %d odd, %d even\n",
+ re->nattempt[0]+re->nattempt[1],re->nattempt[1],re->nattempt[0]);
- fprintf(fplog,"Repl number of exchanges:\n");
- print_ind(fplog,"",re->nrepl,re->ind,NULL);
- print_count(fplog,"",re->nrepl,re->nexchange);
-
- fprintf(fplog,"Repl average number of exchanges:\n");
- for(i=1; i<re->nrepl; i++)
- {
- if (re->nattempt[i%2] == 0)
+ snew(prob,re->nrepl);
+ for(i=1; i<re->nrepl; i++)
{
- prob[i] = 0;
+ if (re->nattempt[i%2] == 0)
+ {
+ prob[i] = 0;
+ }
+ else
+ {
+ prob[i] = re->prob_sum[i]/re->nattempt[i%2];
+ }
}
- else
+ print_ind(fplog,"",re->nrepl,re->ind,NULL);
+ print_prob(fplog,"",re->nrepl,prob);
+
+ fprintf(fplog,"Repl number of exchanges:\n");
+ print_ind(fplog,"",re->nrepl,re->ind,NULL);
+ print_count(fplog,"",re->nrepl,re->nexchange);
+
+ fprintf(fplog,"Repl average number of exchanges:\n");
+ for(i=1; i<re->nrepl; i++)
{
- prob[i] = ((real)re->nexchange[i])/re->nattempt[i%2];
+ if (re->nattempt[i%2] == 0)
+ {
+ prob[i] = 0;
+ }
+ else
+ {
+ prob[i] = ((real)re->nexchange[i])/re->nattempt[i%2];
+ }
}
+ print_ind(fplog,"",re->nrepl,re->ind,NULL);
+ print_prob(fplog,"",re->nrepl,prob);
+ sfree(prob);
+ fprintf(fplog,"\n");
}
- print_ind(fplog,"",re->nrepl,re->ind,NULL);
- print_prob(fplog,"",re->nrepl,prob);
-
- sfree(prob);
-
- fprintf(fplog,"\n");
+ /* print the transition matrix */
+ print_matrix(fplog,"",re->nrepl,re->nmoves,re->nattempt);
}
const gmx_multisim_t *ms,
const t_state *state,
const t_inputrec *ir,
- int nst,int init_seed);
+ int nst, int nmultiex, int init_seed);
/* Should only be called on the master nodes */
extern gmx_bool replica_exchange(FILE *fplog,
const t_commrec *cr,
gmx_repl_ex_t re,
- t_state *state,real *ener,
+ t_state *state,gmx_enerdata_t *enerd,
t_state *state_local,
- int step,real time);
+ gmx_large_int_t step,real time);
/* Attempts replica exchange, should be called on all nodes.
* Returns TRUE if this state has been exchanged.
* When running each replica in parallel,
return min1-min2;
}
- static int icomp(const void *p1, const void *p2)
- {
- atom_id *a1=(atom_id *)p1;
- atom_id *a2=(atom_id *)p2;
-
- return (*a1)-(*a2);
- }
-
int n_flexible_constraints(struct gmx_constr *constr)
{
int nflexcon;
/* Set the constraint lengths for the step at which this configuration
* is meant to be. The invmasses should not be changed.
*/
- lambda += delta_step*ir->delta_lambda;
+ lambda += delta_step*ir->fepvals->delta_lambda;
}
if (vir != NULL)
settle = &idef->il[F_SETTLE];
if (settle->nr > 0)
{
- nsettle = settle->nr/2;
+ nsettle = settle->nr/4;
switch (econq)
{
case econqCoord:
csettle(constr->settled,
nsettle,settle->iatoms,x[0],xprime[0],
- invdt,v[0],vir!=NULL,rmdr,&error,&vetavar);
+ invdt,v?v[0]:NULL,vir!=NULL,rmdr,&error,&vetavar);
inc_nrnb(nrnb,eNR_SETTLE,nsettle);
if (v != NULL)
{
sprintf(buf,
"\nstep " gmx_large_int_pfmt ": Water molecule starting at atom %d can not be "
"settled.\nCheck for bad contacts and/or reduce the timestep if appropriate.\n",
- step,ddglatnr(cr->dd,settle->iatoms[error*2+1]));
+ step,ddglatnr(cr->dd,settle->iatoms[error*4+1]));
if (fplog)
{
fprintf(fplog,"%s",buf);
{
settle = &idef->il[F_SETTLE];
iO = settle->iatoms[1];
- iH = settle->iatoms[1]+1;
+ iH = settle->iatoms[2];
constr->settled =
settle_init(md->massT[iO],md->massT[iH],
md->invmass[iO],md->invmass[iH],
constr_recur(&at2con,molt->ilist,iparams,
TRUE,at,0,1+ir->nProjOrder,path,r0,r1,&r2maxB,&count);
}
- lam0 = ir->init_lambda;
+ lam0 = ir->fepvals->init_lambda;
if (EI_DYNAMICS(ir->eI))
- lam0 += ir->init_step*ir->delta_lambda;
+ lam0 += ir->init_step*ir->fepvals->delta_lambda;
rmax = (1 - lam0)*sqrt(r2maxA) + lam0*sqrt(r2maxB);
if (EI_DYNAMICS(ir->eI)) {
- lam1 = ir->init_lambda + (ir->init_step + ir->nsteps)*ir->delta_lambda;
+ lam1 = ir->fepvals->init_lambda + (ir->init_step + ir->nsteps)*ir->fepvals->delta_lambda;
rmax = max(rmax,(1 - lam1)*sqrt(r2maxA) + lam1*sqrt(r2maxB));
}
}
iloop = gmx_mtop_ilistloop_init(mtop);
while (gmx_mtop_ilistloop_next(iloop,&ilist,&nmol))
{
- for (i=0; i<ilist[F_SETTLE].nr; i+=2)
+ for (i=0; i<ilist[F_SETTLE].nr; i+=4)
{
if (settle_type == -1)
{
return bInterCG;
}
+
+/* helper functions for andersen temperature control, because the
+ * gmx_constr construct is only defined in constr.c. Return the list
+ * of blocks (get_sblock) and the number of blocks (get_nblocks). */
+
+extern int *get_sblock(struct gmx_constr *constr)
+{
+ return constr->sblock;
+}
+
+extern int get_nblocks(struct gmx_constr *constr)
+{
+ return constr->nblocks;
+}
#ifdef GMX_LIB_MPI
#include <mpi.h>
#endif
-#ifdef GMX_THREADS
+#ifdef GMX_THREAD_MPI
#include "tmpi.h"
#endif
#ifdef GMX_MPI
MPI_Comm mpi_comm_mygroup;
#endif
+ int omp_nthreads;
} gmx_wallcycle_t_t;
/* Each name should not exceed 19 characters */
static const char *wcn[ewcNR] =
-{ "Run", "Step", "PP during PME", "Domain decomp.", "DD comm. load", "DD comm. bounds", "Vsite constr.", "Send X to PME", "Comm. coord.", "Neighbor search", "Born radii", "Force", "Wait + Comm. F", "PME mesh", "PME redist. X/F", "PME spread/gather", "PME 3D-FFT", "PME solve", "Wait + Comm. X/F", "Wait + Recv. PME F", "Vsite spread", "Write traj.", "Update", "Constraints", "Comm. energies", "Test" };
+{ "Run", "Step", "PP during PME", "Domain decomp.", "DD comm. load", "DD comm. bounds", "Vsite constr.", "Send X to PME", "Comm. coord.", "Neighbor search", "Born radii", "Force", "Wait + Comm. F", "PME mesh", "PME redist. X/F", "PME spread/gather", "PME 3D-FFT", "PME 3D-FFT Comm.", "PME solve", "Wait + Comm. X/F", "Wait + Recv. PME F", "Vsite spread", "Write traj.", "Update", "Constraints", "Comm. energies", "Enforced rotation", "Add rot. forces", "Test" };
gmx_bool wallcycle_have_counter(void)
{
return gmx_cycles_have_counter();
}
-gmx_wallcycle_t wallcycle_init(FILE *fplog,int resetstep,t_commrec *cr)
+gmx_wallcycle_t wallcycle_init(FILE *fplog,int resetstep,t_commrec *cr, int omp_nthreads)
{
gmx_wallcycle_t wc;
wc->wc_depth = 0;
wc->ewc_prev = -1;
wc->reset_counters = resetstep;
+ wc->omp_nthreads = omp_nthreads;
#ifdef GMX_MPI
if (PAR(cr) && getenv("GMX_CYCLE_BARRIER") != NULL)
snew(wc->wcc,ewcNR);
if (getenv("GMX_CYCLE_ALL") != NULL)
{
-/*#ifndef GMX_THREADS*/
+/*#ifndef GMX_THREAD_MPI*/
if (fplog)
{
fprintf(fplog,"\nWill time all the code during the run\n\n");
return wc;
}
-void wallcycle_destroy(gmx_wallcycle_t wc)
-{
- if (wc == NULL)
- {
- return;
- }
-
- if (wc->wcc != NULL)
- {
- sfree(wc->wcc);
- }
- if (wc->wcc_all != NULL)
- {
- sfree(wc->wcc_all);
- }
- sfree(wc);
-}
-
static void wallcycle_all_start(gmx_wallcycle_t wc,int ewc,gmx_cycles_t cycle)
{
wc->ewc_prev = ewc;
}
}
+static gmx_bool pme_subdivision(int ewc)
+{
+ return (ewc >= ewcPME_REDISTXF && ewc <= ewcPME_SOLVE);
+}
+
void wallcycle_sum(t_commrec *cr, gmx_wallcycle_t wc,double cycles[])
{
wallcc_t *wcc;
wcc = wc->wcc;
+ if (wc->omp_nthreads>1)
+ {
+ for(i=0; i<ewcNR; i++)
+ {
+ if (pme_subdivision(i) || i==ewcPMEMESH || (i==ewcRUN && cr->duty == DUTY_PME))
+ {
+ wcc[i].c *= wc->omp_nthreads;
+ }
+ }
+ }
+
if (wcc[ewcDDCOMMLOAD].n > 0)
{
wcc[ewcDOMDEC].c -= wcc[ewcDDCOMMLOAD].c;
{
wcc[ewcDOMDEC].c -= wcc[ewcDDCOMMBOUND].c;
}
+ if (wcc[ewcPME_FFTCOMM].n > 0)
+ {
+ wcc[ewcPME_FFT].c -= wcc[ewcPME_FFTCOMM].c;
+ }
+
if (cr->npmenodes == 0)
{
/* All nodes do PME (or no PME at all) */
}
}
-static gmx_bool subdivision(int ewc)
-{
- return (ewc >= ewcPME_REDISTXF && ewc <= ewcPME_SOLVE);
-}
void wallcycle_print(FILE *fplog, int nnodes, int npme, double realtime,
gmx_wallcycle_t wc, double cycles[])
npme = nnodes;
}
tot = cycles[ewcRUN];
+ /* PME part has to be multiplied with number of threads */
+ if (npme == 0)
+ {
+ tot += cycles[ewcPMEMESH]*(wc->omp_nthreads-1);
+ }
/* Conversion factor from cycles to seconds */
if (tot > 0)
{
- c2t = nnodes*realtime/tot;
+ c2t = (npp+npme*wc->omp_nthreads)*realtime/tot;
}
else
{
sum = 0;
for(i=ewcPPDURINGPME+1; i<ewcNR; i++)
{
- if (!subdivision(i))
+ if (!pme_subdivision(i))
{
print_cycles(fplog,c2t,wcn[i],
(i==ewcPMEMESH || i==ewcPMEWAITCOMM) ? npme : npp,
fprintf(fplog,"%s\n",myline);
for(i=ewcPPDURINGPME+1; i<ewcNR; i++)
{
- if (subdivision(i))
+ if (pme_subdivision(i))
{
print_cycles(fplog,c2t,wcn[i],
- (i>=ewcPMEMESH || i<=ewcPME_SOLVE) ? npme : npp,
+ (i>=ewcPMEMESH && i<=ewcPME_SOLVE) ? npme : npp,
wc->wcc[i].n,cycles[i],tot);
}
}
#include "smalloc.h"
#include "macros.h"
-#include "math.h"
+#include <math.h>
#include "xvgr.h"
#include "copyrite.h"
#include "statutil.h"
"With options [TT]-lt[tt] and [TT]-dist[tt] the number of contacts",
"of all atoms in group 2 that are closer than a certain distance",
"to the center of mass of group 1 are plotted as a function of the time",
- "that the contact was continuously present.[PAR]",
+ "that the contact was continuously present. The [TT]-intra[tt] switch enables",
+ "calculations of intramolecular distances avoiding distance calculation to its",
+ "periodic images. For a proper function, the molecule in the input trajectory",
+ "should be whole (e.g. by preprocessing with [TT]trjconv -pbc[tt]) or a matching",
+ "topology should be provided. The [TT]-intra[tt] switch will only give",
+ "meaningful results for intramolecular and not intermolecular distances.[PAR]",
"Other programs that calculate distances are [TT]g_mindist[tt]",
"and [TT]g_bond[tt]."
};
atom_id aid;
int ngrps; /* the number of index groups */
- atom_id **index,max; /* the index for the atom numbers */
+ atom_id **index, natoms_groups; /* the index for the atom numbers */
int *isize; /* the size of each group */
char **grpname; /* the name of each group */
rvec *com;
real *mass;
FILE *fp=NULL,*fplt=NULL;
- gmx_bool bCutoff,bPrintDist,bLifeTime;
+ gmx_bool bCutoff,bPrintDist,bLifeTime,bIntra=FALSE;
t_pbc *pbc;
int *contact_time=NULL,*ccount=NULL,ccount_nalloc=0,sum;
char buf[STRLEN];
static real cut=0;
- static t_pargs pa[] = {
+ t_pargs pa[] = {
+ { "-intra", FALSE, etBOOL, {&bIntra},
+ "Calculate distances without considering periodic boundaries, e.g. intramolecular." },
{ "-dist", FALSE, etREAL, {&cut},
"Print all atoms in group 2 closer than dist to the center of mass of group 1" }
};
get_index(&top->atoms,ftp2fn(efNDX,NFILE,fnm),ngrps,isize,index,grpname);
/* calculate mass */
- max=0;
+ natoms_groups=0;
snew(mass,ngrps);
for(g=0;(g<ngrps);g++) {
mass[g]=0;
for(i=0;(i<isize[g]);i++) {
- if (index[g][i]>max)
- max=index[g][i];
+ if (index[g][i]>natoms_groups)
+ natoms_groups=index[g][i];
if (index[g][i] >= top->atoms.nr)
gmx_fatal(FARGS,"Atom number %d, item %d of group %d, is larger than number of atoms in the topolgy (%d)\n",index[g][i]+1,i+1,g+1,top->atoms.nr+1);
mass[g]+=top->atoms.atom[index[g][i]].m;
}
}
+ /* The number is one more than the highest index */
+ natoms_groups++;
natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
t0 = t;
- if (max>=natoms)
- gmx_fatal(FARGS,"Atom number %d in an index group is larger than number of atoms in the trajectory (%d)\n",(int)max+1,natoms);
+ if (natoms_groups>natoms)
+ gmx_fatal(FARGS,"Atom number %d in an index group is larger than number of atoms in the trajectory (%d)\n",(int)natoms_groups,natoms);
if (!bCutoff) {
/* open output file */
else
pbc = NULL;
- gpbc = gmx_rmpbc_init(&top->idef,ePBC,max,box);
+ gpbc = gmx_rmpbc_init(&top->idef,ePBC,natoms,box);
do {
/* initialisation for correct distance calculations */
if (pbc) {
set_pbc(pbc,ePBC,box);
/* make molecules whole again */
- gmx_rmpbc(gpbc,max,box,x);
+ gmx_rmpbc(gpbc,natoms,box,x);
}
/* calculate center of masses */
for(g=0;(g<ngrps);g++) {
/* write to output */
fprintf(fp,"%12.7f ",t);
for(g=0;(g<ngrps/2);g++) {
- if (pbc)
+ if (pbc && (!bIntra))
pbc_dx(pbc,com[2*g],com[2*g+1],dx);
else
rvec_sub(com[2*g],com[2*g+1],dx);
} else {
for(i=0;(i<isize[1]);i++) {
j=index[1][i];
- if (pbc)
+ if (pbc && (!bIntra))
pbc_dx(pbc,x[j],com[0],dx);
else
rvec_sub(x[j],com[0],dx);
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff