Default values are given in parentheses, or listed first among
choices. The first option in the list is always the default
-option. Units are given in square brackets The difference between a
+option. Units are given in square brackets. The difference between a
dash and an underscore is ignored.
A :ref:`sample mdp file <mdp>` is available. This should be
twice as fast as with a Berendsen thermostat with the same
:mdp:`tau-t`.
- .. mdp-value:: sd2
-
- This used to be the default sd integrator, but is now
- deprecated. Four Gaussian random numbers are required per
- coordinate per step. With constraints, the temperature will be
- slightly too high.
-
.. mdp-value:: bd
An Euler integrator for Brownian or position Langevin dynamics,
molecule. When :mdp:`nstlist` is larger than one,
:mdp:`nstlist` insertions are performed in a sphere with radius
:mdp:`rtpi` around a the same random location using the same
- neighborlist (and the same long-range energy when :mdp:`rvdw`
- or :mdp:`rcoulomb` > :mdp:`rlist`, 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
+ neighborlist. Since neighborlist construction is expensive,
+ one can perform several extra insertions with the same list
+ almost for free. The random seed is set with
:mdp:`ld-seed`. The temperature for the Boltzmann weighting is
set with :mdp:`ref-t`, this should match the temperature of the
simulation of the original trajectory. Dispersion correction is
(100)
number of steps that elapse between calculating the energies, 0 is
- never. This option is only relevant with dynamics. With a
- twin-range cut-off setup :mdp:`nstcalcenergy` should be equal to
- or a multiple of :mdp:`nstlist`. This option affects the
+ never. This option is only relevant with dynamics. This option affects the
performance in parallel simulations, because calculating energies
requires global communication between all processes which can
become a bottleneck at high parallelization.
automatically set based on :mdp:`verlet-buffer-tolerance`,
unless this is set to -1, in which case :mdp:`rlist` will be
used. This option has an explicit, exact cut-off at :mdp:`rvdw`
- equal to :mdp:`rcoulomb`. Currently only cut-off,
- reaction-field, PME electrostatics and plain LJ are
- supported. Some :ref:`gmx mdrun` functionality is not yet
+ equal to :mdp:`rcoulomb`, unless PME or Ewald is used, in which
+ case :mdp:`rcoulomb` > :mdp:`rvdw` is allowed. Currently only
+ cut-off, reaction-field, PME or Ewald electrostatics and plain
+ LJ are supported. Some :ref:`gmx mdrun` functionality is not yet
supported with the :mdp:`Verlet` scheme, but :ref:`gmx grompp`
checks for this. Native GPU acceleration is only supported with
:mdp:`Verlet`. With GPU-accelerated PME or with separate PME
.. mdp-value:: >0
- Frequency to update the neighbor list (and the long-range
- forces, when using twin-range cut-offs). When this is 0, the
+ Frequency to update the neighbor list. When this is 0, the
neighbor list is made only once. With energy minimization the
neighborlist will be updated for every energy evaluation when
:mdp:`nstlist` is greater than 0. With :mdp:`Verlet` and
Unused.
-.. mdp:: nstcalclr
-
- (-1) \[steps\]
- Controls the period between calculations of long-range forces when
- using the group cut-off scheme.
-
- .. mdp-value:: 1
-
- Calculate the long-range forces every single step. This is
- useful to have separate neighbor lists with buffers for
- electrostatics and Van der Waals interactions, and in particular
- it makes it possible to have the Van der Waals cutoff longer
- than electrostatics (useful *e.g.* with PME). However, there is
- no point in having identical long-range cutoffs for both
- interaction forms and update them every step - then it will be
- slightly faster to put everything in the short-range list.
-
- .. mdp-value:: >1
-
- Calculate the long-range forces every :mdp:`nstcalclr` steps
- and use a multiple-time-step integrator to combine forces. This
- can now be done more frequently than :mdp:`nstlist` since the
- lists are stored, and it might be a good idea *e.g.* for Van der
- Waals interactions that vary slower than electrostatics.
-
- .. mdp-value:: \-1
-
- Calculate long-range forces on steps where neighbor searching is
- performed. While this is the default value, you might want to
- consider updating the long-range forces more frequently.
-
- Note that twin-range force evaluation might be enabled
- automatically by PP-PME load balancing. This is done in order to
- maintain the chosen Van der Waals interaction radius even if the
- load balancing is changing the electrostatics cutoff. If the
- :ref:`mdp` file already specifies twin-range interactions (*e.g.* to
- evaluate Lennard-Jones interactions with a longer cutoff than
- the PME electrostatics every 2-3 steps), the load balancing will
- have also a small effect on Lennard-Jones, since the short-range
- cutoff (inside which forces are evaluated every step) is
- changed.
-
.. mdp:: ns-type
.. mdp-value:: grid
:mdp:`verlet-buffer-tolerance` option and the value of
:mdp:`rlist` is ignored.
-.. mdp:: rlistlong
-
- (-1) \[nm\]
- Cut-off distance for the long-range neighbor list. This parameter
- is only relevant for a twin-range cut-off setup with switched
- potentials. In that case a buffer region is required to account for
- the size of charge groups. In all other cases this parameter is
- automatically set to the longest cut-off distance.
-
Electrostatics
^^^^^^^^^^^^^^
.. mdp-value:: Cut-off
- Twin range cut-offs with neighborlist cut-off :mdp:`rlist` and
- Coulomb cut-off :mdp:`rcoulomb`, where :mdp:`rcoulomb` >=
- :mdp:`rlist`.
+ Plain cut-off with neighborlist radius :mdp:`rlist` and
+ Coulomb cut-off :mdp:`rcoulomb`, where :mdp:`rlist` >=
+ :mdp:`rcoulomb`.
.. mdp-value:: Ewald
.. mdp-value:: Reaction-Field
Reaction field electrostatics with Coulomb cut-off
- :mdp:`rcoulomb`, where :mdp:`rcoulomb` >= :mdp:`rlist`. The
+ :mdp:`rcoulomb`, where :mdp:`rlist` >= :mdp:`rvdw`. The
dielectric constant beyond the cut-off is
:mdp:`epsilon-rf`. The dielectric constant can be set to
infinity by setting :mdp:`epsilon-rf` =0.
.. mdp-value:: Generalized-Reaction-Field
Generalized reaction field with Coulomb cut-off
- :mdp:`rcoulomb`, where :mdp:`rcoulomb` >= :mdp:`rlist`. The
+ :mdp:`rcoulomb`, where :mdp:`rlist` >= :mdp:`rcoulomb`. The
dielectric constant beyond the cut-off is
:mdp:`epsilon-rf`. The ionic strength is computed from the
number of charged (*i.e.* with non zero charge) charge
:mdp:`Reaction-Field-zero` computationally more expensive than
normal reaction-field.
- .. mdp-value:: Reaction-Field-nec
-
- The same as :mdp-value:`coulombtype=Reaction-Field`, but
- implemented as in |Gromacs| versions before 3.3. No
- reaction-field correction is applied to excluded atom pairs and
- self pairs. The 1-4 interactions are calculated using a
- reaction-field. The missing correction due to the excluded pairs
- that do not have a 1-4 interaction is up to a few percent of the
- total electrostatic energy and causes a minor difference in the
- forces and the pressure.
-
.. mdp-value:: Shift
Analogous to :mdp-value:`vdwtype=Shift` for :mdp:`vdwtype`. You
.. mdp:: pcoupltype
+ Specifies the kind of isotropy of the pressure coupling used. Each
+ kind takes one or more values for :mdp:`compressibility` and
+ :mdp:`ref-p`. Only a single value is permitted for :mdp:`tau-p`.
+
.. mdp-value:: isotropic
Isotropic pressure coupling with time constant
- :mdp:`tau-p`. The compressibility and reference pressure are
- set with :mdp:`compressibility` and :mdp:`ref-p`, one value is
- needed.
+ :mdp:`tau-p`. One value each for :mdp:`compressibility` and
+ :mdp:`ref-p` is required.
.. mdp-value:: semiisotropic
Pressure coupling which is isotropic in the ``x`` and ``y``
direction, but different in the ``z`` direction. This can be
- useful for membrane simulations. 2 values are needed for ``x/y``
- and ``z`` directions respectively.
+ useful for membrane simulations. Two values each for
+ :mdp:`compressibility` and :mdp:`ref-p` are required, for
+ ``x/y`` and ``z`` directions respectively.
.. mdp-value:: anisotropic
.. mdp:: tau-p
(1) \[ps\]
- time constant for coupling
+ The time constant for pressure coupling (one value for all
+ directions).
.. mdp:: compressibility
\[bar^-1\]
- compressibility (NOTE: this is now really in bar-1) For water at 1
- atm and 300 K the compressibility is 4.5e-5 bar^-1.
+ The compressibility (NOTE: this is now really in bar^-1) For water at 1
+ atm and 300 K the compressibility is 4.5e-5 bar^-1. The number of
+ required values is implied by :mdp:`pcoupltype`.
.. mdp:: ref-p
\[bar\]
- reference pressure for coupling
+ The reference pressure for coupling. The number of required values
+ is implied by :mdp:`pcoupltype`.
.. mdp:: refcoord-scaling
.. mdp:: lincs-warnangle
- (30) \[degrees\]
+ (30) \[deg\]
maximum angle that a bond can rotate before LINCS will complain
.. mdp:: morse
(1e-6)
the relative constraint tolerance for constraint pulling
-.. mdp:: pull-print-com1
+.. mdp:: pull-print-com
.. mdp-value:: no
- do not print the COM of the first group in each pull coordinate
+ do not print the COM for any group
.. mdp-value:: yes
- print the COM of the first group in each pull coordinate
-
-.. mdp:: pull-print-com2
-
- .. mdp-value:: no
-
- do not print the COM of the second group in each pull coordinate
-
- .. mdp-value:: yes
-
- print the COM of the second group in each pull coordinate
+ print the COM of all groups for all pull coordinates
.. mdp:: pull-print-ref-value
.. mdp-value:: flat-bottom
- At distances beyond :mdp:`pull-coord1-init` a harmonic potential
+ At distances above :mdp:`pull-coord1-init` a harmonic potential
+ is applied, otherwise no potential is applied.
+
+ .. mdp-value:: flat-bottom-high
+
+ At distances below :mdp:`pull-coord1-init` a harmonic potential
is applied, otherwise no potential is applied.
+ .. mdp-value:: external-potential
+
+ An external potential that needs to be provided by another
+ module.
+
+.. mdp:: pull-coord1-potential-provider
+
+ The name of the external module that provides the potential for
+ the case where :mdp:`pull-coord1-type` is external-potential.
+
.. mdp:: pull-coord1-geometry
.. mdp-value:: distance
component. This geometry is not supported with constraint
pulling.
+ .. mdp-value:: angle
+
+ Pull along an angle defined by four groups. The angle is
+ defined as the angle between two vectors: the vector connecting
+ the COM of the first group to the COM of the second group and
+ the vector connecting the COM of the third group to the COM of
+ the fourth group.
+
+ .. mdp-value:: angle-axis
+
+ As :mdp-value:`angle` but the second vector is given by :mdp:`pull-coord1-vec`.
+ Thus, only the two groups that define the first vector need to be given.
+
+ .. mdp-value:: dihedral
+
+ Pull along a dihedral angle defined by six groups. These pairwise
+ define three vectors: the vector connecting the COM of group 1
+ to the COM of group 2, the COM of group 3 to the COM of group 4,
+ and the COM of group 5 to the COM group 6. The dihedral angle is
+ then defined as the angle between two planes: the plane spanned by the
+ the two first vectors and the plane spanned the two last vectors.
+
+
.. mdp:: pull-coord1-groups
- The two groups indices should be given on which this pull
- coordinate will operate. The first index can be 0, in which case an
+ The group indices on which this pull coordinate will operate.
+ The number of group indices required is geometry dependent.
+ The first index can be 0, in which case an
absolute reference of :mdp:`pull-coord1-origin` is used. With an
absolute reference the system is no longer translation invariant
and one should think about what to do with the center of mass
- motion. Note that (only) for :mdp:`pull-coord1-geometry` =
- :mdp-value:`direction-relative` four groups are required.
+ motion.
.. mdp:: pull-coord1-dim
.. mdp:: pull-coord1-init
- (0.0) \[nm\]
+ (0.0) \[nm\] / \[deg\]
The reference distance at t=0.
.. mdp:: pull-coord1-rate
- (0) \[nm/ps\]
+ (0) \[nm/ps\] / \[deg/ps\]
The rate of change of the reference position.
.. mdp:: pull-coord1-k
- (0) \[kJ mol-1 nm-2\] / \[kJ mol-1 nm-1\]
+ (0) \[kJ mol-1 nm-2\] / \[kJ mol-1 nm-1\] / \[kJ mol-1 rad-2\] / \[kJ mol-1 rad-1\]
The force constant. For umbrella pulling this is the harmonic force
- constant in kJ mol-1 nm-2. For constant force pulling this is the
+ constant in kJ mol-1 nm-2 (or kJ mol-1 rad-2 for angles). For constant force pulling this is the
force constant of the linear potential, and thus the negative (!)
- of the constant force in kJ mol-1 nm-1.
+ of the constant force in kJ mol-1 nm-1 (or kJ mol-1 rad-1 for angles).
+ Note that for angles the force constant is expressed in terms of radians
+ (while :mdp:`pull-coord1-init` and :mdp:`pull-coord1-rate` are expressed in degrees).
.. mdp:: pull-coord1-kB
- (pull-k1) \[kJ mol-1 nm-2\] / \[kJ mol-1 nm-1\]
+ (pull-k1) \[kJ mol-1 nm-2\] / \[kJ mol-1 nm-1\] / \[kJ mol-1 rad-2\] / \[kJ mol-1 rad-1\]
As :mdp:`pull-coord1-k`, but for state B. This is only used when
:mdp:`free-energy` is turned on. The force constant is then (1 -
lambda) * :mdp:`pull-coord1-k` + lambda * :mdp:`pull-coord1-kB`.
calculations are done.
-Adaptive Resolution Simulation
-^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Computational Electrophysiology
+^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
+Use these options to switch on and control ion/water position exchanges in "Computational
+Electrophysiology" simulation setups. (See the `reference manual`_ for details).
-.. mdp:: adress
+.. mdp:: swapcoords
- (no)
- Decide whether the AdResS feature is turned on.
+ .. mdp-value:: no
-.. mdp:: adress-type
+ Do not enable ion/water position exchanges.
- .. mdp-value:: Off
+ .. mdp-value:: X ; Y ; Z
- 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.
+ Allow for ion/water position exchanges along the chosen direction.
+ In a typical setup with the membranes parallel to the x-y plane,
+ ion/water pairs need to be exchanged in Z direction to sustain the
+ requested ion concentrations in the compartments.
- .. mdp-value:: Constant
+.. mdp:: swap-frequency
- Do an AdResS simulation with a constant weight,
- :mdp:`adress-const-wf` defines the value of the weight
+ (1) The swap attempt frequency, i.e. every how many time steps the ion counts
+ per compartment are determined and exchanges made if necessary.
+ Normally it is not necessary to check at every time step.
+ For typical Computational Electrophysiology setups, a value of about 100 is
+ sufficient and yields a negligible performance impact.
- .. mdp-value:: XSplit
+.. mdp:: split-group0
- 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.
+ Name of the index group of the membrane-embedded part of channel #0.
+ The center of mass of these atoms defines one of the compartment boundaries
+ and should be chosen such that it is near the center of the membrane.
- .. mdp-value:: Sphere
+.. mdp:: split-group1
- Do an AdResS simulation with spherical explicit zone.
+ Channel #1 defines the position of the other compartment boundary.
-.. mdp:: adress-const-wf
+.. mdp:: massw-split0
- (1)
- Provides the weight for a constant weight simulation
- (:mdp:`adress-type` =Constant)
+ (no) Defines whether or not mass-weighting is used to calculate the split group center.
-.. mdp:: adress-ex-width
+ .. mdp-value:: no
- (0)
- Width of the explicit zone, measured from
- :mdp:`adress-reference-coords`.
+ Use the geometrical center.
-.. mdp:: adress-hy-width
+ .. mdp-value:: yes
- (0)
- Width of the hybrid zone.
+ Use the center of mass.
-.. mdp:: adress-reference-coords
+.. mdp:: massw-split1
- (0,0,0)
- Position of the center of the explicit zone. Periodic boundary
- conditions apply for measuring the distance from it.
+ (no) As above, but for split-group #1.
-.. mdp:: adress-cg-grp-names
+.. mdp:: solvent-group
- 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.
+ Name of the index group of solvent molecules.
-.. mdp:: adress-site
+.. mdp:: coupl-steps
- The mapping point from which the weight is calculated.
+ (\10) Average the number of ions per compartment over these many swap attempt steps.
+ This can be used to prevent that ions near a compartment boundary
+ (diffusing through a channel, e.g.) lead to unwanted back and forth swaps.
- .. mdp-value:: COM
+.. mdp:: iontypes
- The weight is calculated from the center of mass of each charge group.
+ (1) The number of different ion types to be controlled. These are during the
+ simulation exchanged with solvent molecules to reach the desired reference numbers.
- .. mdp-value:: COG
+.. mdp:: iontype0-name
- The weight is calculated from the center of geometry of each charge group.
+ Name of the first ion type.
- .. mdp-value:: Atom
+.. mdp:: iontype0-in-A
- The weight is calculated from the position of 1st atom of each charge group.
+ (-1) Requested (=reference) number of ions of type 0 in compartment A.
+ The default value of -1 means: use the number of ions as found in time step 0
+ as reference value.
- .. mdp-value:: AtomPerAtom
+.. mdp:: iontype0-in-B
- The weight is calculated from the position of each individual atom.
+ (-1) Reference number of ions of type 0 for compartment B.
-.. mdp:: adress-interface-correction
+.. mdp:: bulk-offsetA
- .. mdp-value:: Off
+ (0.0) Offset of the first swap layer from the compartment A midplane.
+ By default (i.e. bulk offset = 0.0), ion/water exchanges happen between layers
+ at maximum distance (= bulk concentration) to the split group layers. However,
+ an offset b (-1.0 < b < +1.0) can be specified to offset the bulk layer from the middle at 0.0
+ towards one of the compartment-partitioning layers (at +/- 1.0).
- Do not apply any interface correction.
+.. mdp:: bulk-offsetB
- .. mdp-value:: thermoforce
+ (0.0) Offset of the other swap layer from the compartment B midplane.
- Apply thermodynamic force interface correction. The table can be
- specified using the ``-tabletf`` option of :ref:`gmx mdrun`. The
- table should contain the potential and force (acting on
- molecules) as function of the distance from
- :mdp:`adress-reference-coords`.
-.. mdp:: adress-tf-grp-names
+.. mdp:: threshold
- The names of the energy groups to which the thermoforce is applied
- if enabled in :mdp:`adress-interface-correction`. If no group is
- given the default table is applied.
+ (\1) Only swap ions if threshold difference to requested count is reached.
-.. mdp:: adress-ex-forcecap
+.. mdp:: cyl0-r
- (0)
- Cap the force in the hybrid region, useful for big molecules. 0
- disables force capping.
+ (2.0) \[nm\] Radius of the split cylinder #0.
+ Two split cylinders (mimicking the channel pores) can optionally be defined
+ relative to the center of the split group. With the help of these cylinders
+ it can be counted which ions have passed which channel. The split cylinder
+ definition has no impact on whether or not ion/water swaps are done.
+
+.. mdp:: cyl0-up
+
+ (1.0) \[nm\] Upper extension of the split cylinder #0.
+
+.. mdp:: cyl0-down
+
+ (1.0) \[nm\] Lower extension of the split cylinder #0.
+
+.. mdp:: cyl1-r
+
+ (2.0) \[nm\] Radius of the split cylinder #1.
+
+.. mdp:: cyl1-up
+
+ (1.0) \[nm\] Upper extension of the split cylinder #1.
+
+.. mdp:: cyl1-down
+
+ (1.0) \[nm\] Lower extension of the split cylinder #1.
User defined thingies
These you can use if you modify code. You can pass integers and
reals and groups to your subroutine. Check the inputrec definition
- in ``src/gromacs/legacyheaders/types/inputrec.h``
+ in ``src/gromacs/mdtypes/inputrec.h``
+
+Removed features
+^^^^^^^^^^^^^^^^
+
+This feature has been removed from |Gromacs|, but so that old
+:ref:`mdp` and :ref:`tpr` files cannot be mistakenly misused, we still
+parse this option. :ref:`gmx grompp` and :ref:`gmx mdrun` will issue a
+fatal error if this is set.
+
+.. mdp:: adress
+
+ (no)
.. _reference manual: gmx-manual-parent-dir_
--- /dev/null
- * Copyright (c) 2013,2014,2015,2016, by the GROMACS development team, led by
+/*
+ * This file is part of the GROMACS molecular simulation package.
+ *
+ * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
+ * Copyright (c) 2001-2004, The GROMACS development team.
- if (z_ave < 0)
++ * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
+ * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
+ * and including many others, as listed in the AUTHORS file in the
+ * top-level source directory and at http://www.gromacs.org.
+ *
+ * GROMACS is free software; you can redistribute it and/or
+ * modify it under the terms of the GNU Lesser General Public License
+ * as published by the Free Software Foundation; either version 2.1
+ * of the License, or (at your option) any later version.
+ *
+ * GROMACS is distributed in the hope that it will be useful,
+ * but WITHOUT ANY WARRANTY; without even the implied warranty of
+ * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
+ * Lesser General Public License for more details.
+ *
+ * You should have received a copy of the GNU Lesser General Public
+ * License along with GROMACS; if not, see
+ * http://www.gnu.org/licenses, or write to the Free Software Foundation,
+ * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
+ *
+ * If you want to redistribute modifications to GROMACS, please
+ * consider that scientific software is very special. Version
+ * control is crucial - bugs must be traceable. We will be happy to
+ * consider code for inclusion in the official distribution, but
+ * derived work must not be called official GROMACS. Details are found
+ * in the README & COPYING files - if they are missing, get the
+ * official version at http://www.gromacs.org.
+ *
+ * To help us fund GROMACS development, we humbly ask that you cite
+ * the research papers on the package. Check out http://www.gromacs.org.
+ */
+#include "gmxpre.h"
+
+#include <cmath>
+#include <cstring>
+
+#include <algorithm>
+
+#include "gromacs/commandline/pargs.h"
+#include "gromacs/commandline/viewit.h"
+#include "gromacs/fileio/confio.h"
+#include "gromacs/fileio/trxio.h"
+#include "gromacs/fileio/xvgr.h"
+#include "gromacs/gmxana/cmat.h"
+#include "gromacs/gmxana/gmx_ana.h"
+#include "gromacs/gmxana/gstat.h"
+#include "gromacs/math/functions.h"
+#include "gromacs/math/utilities.h"
+#include "gromacs/math/vec.h"
+#include "gromacs/pbcutil/pbc.h"
+#include "gromacs/pbcutil/rmpbc.h"
+#include "gromacs/topology/index.h"
+#include "gromacs/topology/topology.h"
+#include "gromacs/trajectory/trajectoryframe.h"
+#include "gromacs/utility/arraysize.h"
+#include "gromacs/utility/cstringutil.h"
+#include "gromacs/utility/fatalerror.h"
+#include "gromacs/utility/futil.h"
+#include "gromacs/utility/gmxassert.h"
+#include "gromacs/utility/smalloc.h"
+
+/****************************************************************************/
+/* This program calculates the order parameter per atom for an interface or */
+/* bilayer, averaged over time. */
+/* S = 1/2 * (3 * cos(i)cos(j) - delta(ij)) */
+/* It is assumed that the order parameter with respect to a box-axis */
+/* is calculated. In that case i = j = axis, and delta(ij) = 1. */
+/* */
+/* Peter Tieleman, April 1995 */
+/* P.J. van Maaren, November 2005 Added tetrahedral stuff */
+/****************************************************************************/
+
+static void find_nearest_neighbours(int ePBC,
+ int natoms, matrix box,
+ rvec x[], int maxidx, int index[],
+ real *sgmean, real *skmean,
+ int nslice, int slice_dim,
+ real sgslice[], real skslice[],
+ gmx_rmpbc_t gpbc)
+{
+ int ix, jx, nsgbin, *sgbin;
+ int i, ibin, j, k, l, n, *nn[4];
+ rvec dx, rj, rk, urk, urj;
+ real cost, cost2, *sgmol, *skmol, rmean, rmean2, r2, box2, *r_nn[4];
+ t_pbc pbc;
+ int sl_index;
+ int *sl_count;
+ real onethird = 1.0/3.0;
+ /* dmat = init_mat(maxidx, FALSE); */
+ box2 = box[XX][XX] * box[XX][XX];
+ snew(sl_count, nslice);
+ for (i = 0; (i < 4); i++)
+ {
+ snew(r_nn[i], natoms);
+ snew(nn[i], natoms);
+
+ for (j = 0; (j < natoms); j++)
+ {
+ r_nn[i][j] = box2;
+ }
+ }
+
+ snew(sgmol, maxidx);
+ snew(skmol, maxidx);
+
+ /* Must init pbc every step because of pressure coupling */
+ set_pbc(&pbc, ePBC, box);
+
+ gmx_rmpbc(gpbc, natoms, box, x);
+
+ nsgbin = 2001; // Calculated as (1 + 1/0.0005)
+ snew(sgbin, nsgbin);
+
+ *sgmean = 0.0;
+ *skmean = 0.0;
+ l = 0;
+ for (i = 0; (i < maxidx); i++) /* loop over index file */
+ {
+ ix = index[i];
+ for (j = 0; (j < maxidx); j++)
+ {
+ if (i == j)
+ {
+ continue;
+ }
+
+ jx = index[j];
+
+ pbc_dx(&pbc, x[ix], x[jx], dx);
+ r2 = iprod(dx, dx);
+
+ /* set_mat_entry(dmat,i,j,r2); */
+
+ /* determine the nearest neighbours */
+ if (r2 < r_nn[0][i])
+ {
+ r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
+ r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
+ r_nn[1][i] = r_nn[0][i]; nn[1][i] = nn[0][i];
+ r_nn[0][i] = r2; nn[0][i] = j;
+ }
+ else if (r2 < r_nn[1][i])
+ {
+ r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
+ r_nn[2][i] = r_nn[1][i]; nn[2][i] = nn[1][i];
+ r_nn[1][i] = r2; nn[1][i] = j;
+ }
+ else if (r2 < r_nn[2][i])
+ {
+ r_nn[3][i] = r_nn[2][i]; nn[3][i] = nn[2][i];
+ r_nn[2][i] = r2; nn[2][i] = j;
+ }
+ else if (r2 < r_nn[3][i])
+ {
+ r_nn[3][i] = r2; nn[3][i] = j;
+ }
+ }
+
+
+ /* calculate mean distance between nearest neighbours */
+ rmean = 0;
+ for (j = 0; (j < 4); j++)
+ {
+ r_nn[j][i] = std::sqrt(r_nn[j][i]);
+ rmean += r_nn[j][i];
+ }
+ rmean /= 4;
+
+ n = 0;
+ sgmol[i] = 0.0;
+ skmol[i] = 0.0;
+
+ /* Chau1998a eqn 3 */
+ /* angular part tetrahedrality order parameter per atom */
+ for (j = 0; (j < 3); j++)
+ {
+ for (k = j+1; (k < 4); k++)
+ {
+ pbc_dx(&pbc, x[ix], x[index[nn[k][i]]], rk);
+ pbc_dx(&pbc, x[ix], x[index[nn[j][i]]], rj);
+
+ unitv(rk, urk);
+ unitv(rj, urj);
+
+ cost = iprod(urk, urj) + onethird;
+ cost2 = cost * cost;
+
+ /* sgmol[i] += 3*cost2/32; */
+ sgmol[i] += cost2;
+
+ /* determine distribution */
+ ibin = static_cast<int>(nsgbin * cost2);
+ if (ibin < nsgbin)
+ {
+ sgbin[ibin]++;
+ }
+ /* printf("%d %d %f %d %d\n", j, k, cost * cost, ibin, sgbin[ibin]);*/
+ l++;
+ n++;
+ }
+ }
+
+ /* normalize sgmol between 0.0 and 1.0 */
+ sgmol[i] = 3*sgmol[i]/32;
+ *sgmean += sgmol[i];
+
+ /* distance part tetrahedrality order parameter per atom */
+ rmean2 = 4 * 3 * rmean * rmean;
+ for (j = 0; (j < 4); j++)
+ {
+ skmol[i] += (rmean - r_nn[j][i]) * (rmean - r_nn[j][i]) / rmean2;
+ /* printf("%d %f (%f %f %f %f) \n",
+ i, skmol[i], rmean, rmean2, r_nn[j][i], (rmean - r_nn[j][i]) );
+ */
+ }
+
+ *skmean += skmol[i];
+
+ /* Compute sliced stuff */
+ sl_index = static_cast<int>(std::round((1+x[i][slice_dim]/box[slice_dim][slice_dim])*nslice)) % nslice;
+ sgslice[sl_index] += sgmol[i];
+ skslice[sl_index] += skmol[i];
+ sl_count[sl_index]++;
+ } /* loop over entries in index file */
+
+ *sgmean /= maxidx;
+ *skmean /= maxidx;
+
+ for (i = 0; (i < nslice); i++)
+ {
+ if (sl_count[i] > 0)
+ {
+ sgslice[i] /= sl_count[i];
+ skslice[i] /= sl_count[i];
+ }
+ }
+ sfree(sl_count);
+ sfree(sgbin);
+ sfree(sgmol);
+ sfree(skmol);
+ for (i = 0; (i < 4); i++)
+ {
+ sfree(r_nn[i]);
+ sfree(nn[i]);
+ }
+}
+
+
+static void calc_tetra_order_parm(const char *fnNDX, const char *fnTPS,
+ const char *fnTRX, const char *sgfn,
+ const char *skfn,
+ int nslice, int slice_dim,
+ const char *sgslfn, const char *skslfn,
+ const gmx_output_env_t *oenv)
+{
+ FILE *fpsg = NULL, *fpsk = NULL;
+ t_topology top;
+ int ePBC;
+ t_trxstatus *status;
+ int natoms;
+ real t;
+ rvec *xtop, *x;
+ matrix box;
+ real sg, sk;
+ int **index;
+ char **grpname;
+ int i, *isize, ng, nframes;
+ real *sg_slice, *sg_slice_tot, *sk_slice, *sk_slice_tot;
+ gmx_rmpbc_t gpbc = NULL;
+
+
+ read_tps_conf(fnTPS, &top, &ePBC, &xtop, NULL, box, FALSE);
+
+ snew(sg_slice, nslice);
+ snew(sk_slice, nslice);
+ snew(sg_slice_tot, nslice);
+ snew(sk_slice_tot, nslice);
+ ng = 1;
+ /* get index groups */
+ printf("Select the group that contains the atoms you want to use for the tetrahedrality order parameter calculation:\n");
+ snew(grpname, ng);
+ snew(index, ng);
+ snew(isize, ng);
+ get_index(&top.atoms, fnNDX, ng, isize, index, grpname);
+
+ /* Analyze trajectory */
+ natoms = read_first_x(oenv, &status, fnTRX, &t, &x, box);
+ if (natoms > top.atoms.nr)
+ {
+ gmx_fatal(FARGS, "Topology (%d atoms) does not match trajectory (%d atoms)",
+ top.atoms.nr, natoms);
+ }
+ check_index(NULL, ng, index[0], NULL, natoms);
+
+ fpsg = xvgropen(sgfn, "S\\sg\\N Angle Order Parameter", "Time (ps)", "S\\sg\\N",
+ oenv);
+ fpsk = xvgropen(skfn, "S\\sk\\N Distance Order Parameter", "Time (ps)", "S\\sk\\N",
+ oenv);
+
+ /* loop over frames */
+ gpbc = gmx_rmpbc_init(&top.idef, ePBC, natoms);
+ nframes = 0;
+ do
+ {
+ find_nearest_neighbours(ePBC, natoms, box, x, isize[0], index[0],
+ &sg, &sk, nslice, slice_dim, sg_slice, sk_slice, gpbc);
+ for (i = 0; (i < nslice); i++)
+ {
+ sg_slice_tot[i] += sg_slice[i];
+ sk_slice_tot[i] += sk_slice[i];
+ }
+ fprintf(fpsg, "%f %f\n", t, sg);
+ fprintf(fpsk, "%f %f\n", t, sk);
+ nframes++;
+ }
+ while (read_next_x(oenv, status, &t, x, box));
+ close_trj(status);
+ gmx_rmpbc_done(gpbc);
+
+ sfree(grpname);
+ sfree(index);
+ sfree(isize);
+
+ xvgrclose(fpsg);
+ xvgrclose(fpsk);
+
+ fpsg = xvgropen(sgslfn,
+ "S\\sg\\N Angle Order Parameter / Slab", "(nm)", "S\\sg\\N",
+ oenv);
+ fpsk = xvgropen(skslfn,
+ "S\\sk\\N Distance Order Parameter / Slab", "(nm)", "S\\sk\\N",
+ oenv);
+ for (i = 0; (i < nslice); i++)
+ {
+ fprintf(fpsg, "%10g %10g\n", (i+0.5)*box[slice_dim][slice_dim]/nslice,
+ sg_slice_tot[i]/nframes);
+ fprintf(fpsk, "%10g %10g\n", (i+0.5)*box[slice_dim][slice_dim]/nslice,
+ sk_slice_tot[i]/nframes);
+ }
+ xvgrclose(fpsg);
+ xvgrclose(fpsk);
+}
+
+
+/* Print name of first atom in all groups in index file */
+static void print_types(int index[], int a[], int ngrps,
+ char *groups[], const t_topology *top)
+{
+ int i;
+
+ fprintf(stderr, "Using following groups: \n");
+ for (i = 0; i < ngrps; i++)
+ {
+ fprintf(stderr, "Groupname: %s First atomname: %s First atomnr %d\n",
+ groups[i], *(top->atoms.atomname[a[index[i]]]), a[index[i]]);
+ }
+ fprintf(stderr, "\n");
+}
+
+static void check_length(real length, int a, int b)
+{
+ if (length > 0.3)
+ {
+ fprintf(stderr, "WARNING: distance between atoms %d and "
+ "%d > 0.3 nm (%f). Index file might be corrupt.\n",
+ a, b, length);
+ }
+}
+
+void calc_order(const char *fn, int *index, int *a, rvec **order,
+ real ***slOrder, real *slWidth, int nslices, gmx_bool bSliced,
+ gmx_bool bUnsat, const t_topology *top, int ePBC, int ngrps, int axis,
+ gmx_bool permolecule, gmx_bool radial, gmx_bool distcalc, const char *radfn,
+ real ***distvals,
+ const gmx_output_env_t *oenv)
+{
+ /* if permolecule = TRUE, order parameters will be calculed per molecule
+ * and stored in slOrder with #slices = # molecules */
+ rvec *x0, /* coordinates with pbc */
+ *x1, /* coordinates without pbc */
+ dist; /* vector between two atoms */
+ matrix box; /* box (3x3) */
+ t_trxstatus *status;
+ rvec cossum, /* sum of vector angles for three axes */
+ Sx, Sy, Sz, /* the three molecular axes */
+ tmp1, tmp2, /* temp. rvecs for calculating dot products */
+ frameorder; /* order parameters for one frame */
+ real *slFrameorder; /* order parameter for one frame, per slice */
+ real length, /* total distance between two atoms */
+ t, /* time from trajectory */
+ z_ave, z1, z2; /* average z, used to det. which slice atom is in */
+ int natoms, /* nr. atoms in trj */
+ nr_tails, /* nr tails, to check if index file is correct */
+ size = 0, /* nr. of atoms in group. same as nr_tails */
+ i, j, m, k, teller = 0,
+ slice, /* current slice number */
+ nr_frames = 0;
+ int *slCount; /* nr. of atoms in one slice */
+ real sdbangle = 0; /* sum of these angles */
+ gmx_bool use_unitvector = FALSE; /* use a specified unit vector instead of axis to specify unit normal*/
+ rvec direction, com, dref, dvec;
+ int comsize, distsize;
+ int *comidx = NULL, *distidx = NULL;
+ char *grpname = NULL;
+ t_pbc pbc;
+ real arcdist, tmpdist;
+ gmx_rmpbc_t gpbc = NULL;
+
+ /* PBC added for center-of-mass vector*/
+ /* Initiate the pbc structure */
+ std::memset(&pbc, 0, sizeof(pbc));
+
+ if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
+ {
+ gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
+ }
+
+ nr_tails = index[1] - index[0];
+ fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
+ /* take first group as standard. Not rocksolid, but might catch error in index*/
+
+ if (permolecule)
+ {
+ nslices = nr_tails;
+ bSliced = FALSE; /*force slices off */
+ fprintf(stderr, "Calculating order parameters for each of %d molecules\n",
+ nslices);
+ }
+
+ if (radial)
+ {
+ use_unitvector = TRUE;
+ fprintf(stderr, "Select an index group to calculate the radial membrane normal\n");
+ get_index(&top->atoms, radfn, 1, &comsize, &comidx, &grpname);
+ }
+ if (distcalc)
+ {
+ if (grpname != NULL)
+ {
+ sfree(grpname);
+ }
+ fprintf(stderr, "Select an index group to use as distance reference\n");
+ get_index(&top->atoms, radfn, 1, &distsize, &distidx, &grpname);
+ bSliced = FALSE; /*force slices off*/
+ }
+
+ if (use_unitvector && bSliced)
+ {
+ fprintf(stderr, "Warning: slicing and specified unit vectors are not currently compatible\n");
+ }
+
+ snew(slCount, nslices);
+ snew(*slOrder, nslices);
+ for (i = 0; i < nslices; i++)
+ {
+ snew((*slOrder)[i], ngrps);
+ }
+ if (distcalc)
+ {
+ snew(*distvals, nslices);
+ for (i = 0; i < nslices; i++)
+ {
+ snew((*distvals)[i], ngrps);
+ }
+ }
+ snew(*order, ngrps);
+ snew(slFrameorder, nslices);
+ snew(x1, natoms);
+
+ if (bSliced)
+ {
+ *slWidth = box[axis][axis]/nslices;
+ fprintf(stderr, "Box divided in %d slices. Initial width of slice: %f\n",
+ nslices, *slWidth);
+ }
+
+
+#if 0
+ nr_tails = index[1] - index[0];
+ fprintf(stderr, "Number of elements in first group: %d\n", nr_tails);
+ /* take first group as standard. Not rocksolid, but might catch error
+ in index*/
+#endif
+
+ teller = 0;
+
+ gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
+ /*********** Start processing trajectory ***********/
+ do
+ {
+ if (bSliced)
+ {
+ *slWidth = box[axis][axis]/nslices;
+ }
+ teller++;
+
+ set_pbc(&pbc, ePBC, box);
+ gmx_rmpbc_copy(gpbc, natoms, box, x0, x1);
+
+ /* Now loop over all groups. There are ngrps groups, the order parameter can
+ be calculated for grp 1 to grp ngrps - 1. For each group, loop over all
+ atoms in group, which is index[i] to (index[i+1] - 1) See block.h. Of
+ course, in this case index[i+1] -index[i] has to be the same for all
+ groups, namely the number of tails. i just runs over all atoms in a tail,
+ so for DPPC ngrps = 16 and i runs from 1 to 14, including 14
+ */
+
+
+ if (radial)
+ {
+ /*center-of-mass determination*/
+ com[XX] = 0.0; com[YY] = 0.0; com[ZZ] = 0.0;
+ for (j = 0; j < comsize; j++)
+ {
+ rvec_inc(com, x1[comidx[j]]);
+ }
+ svmul(1.0/comsize, com, com);
+ }
+ if (distcalc)
+ {
+ dref[XX] = 0.0; dref[YY] = 0.0; dref[ZZ] = 0.0;
+ for (j = 0; j < distsize; j++)
+ {
+ rvec_inc(dist, x1[distidx[j]]);
+ }
+ svmul(1.0/distsize, dref, dref);
+ if (radial)
+ {
+ pbc_dx(&pbc, dref, com, dvec);
+ unitv(dvec, dvec);
+ }
+ }
+
+ for (i = 1; i < ngrps - 1; i++)
+ {
+ clear_rvec(frameorder);
+
+ size = index[i+1] - index[i];
+ if (size != nr_tails)
+ {
+ gmx_fatal(FARGS, "grp %d does not have same number of"
+ " elements as grp 1\n", i);
+ }
+
+ for (j = 0; j < size; j++)
+ {
+ if (radial)
+ /*create unit vector*/
+ {
+ pbc_dx(&pbc, x1[a[index[i]+j]], com, direction);
+ unitv(direction, direction);
+ /*DEBUG*/
+ /*if (j==0)
+ fprintf(stderr,"X %f %f %f\tcom %f %f %f\tdirection %f %f %f\n",x1[a[index[i]+j]][0],x1[a[index[i]+j]][1],x1[a[index[i]+j]][2],com[0],com[1],com[2],
+ direction[0],direction[1],direction[2]);*/
+ }
+
+ if (bUnsat)
+ {
+ /* Using convention for unsaturated carbons */
+ /* first get Sz, the vector from Cn to Cn+1 */
+ rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], dist);
+ length = norm(dist);
+ check_length(length, a[index[i]+j], a[index[i+1]+j]);
+ svmul(1.0/length, dist, Sz);
+
+ /* this is actually the cosine of the angle between the double bond
+ and axis, because Sz is normalized and the two other components of
+ the axis on the bilayer are zero */
+ if (use_unitvector)
+ {
+ sdbangle += gmx_angle(direction, Sz); /*this can probably be optimized*/
+ }
+ else
+ {
+ sdbangle += std::acos(Sz[axis]);
+ }
+ }
+ else
+ {
+ /* get vector dist(Cn-1,Cn+1) for tail atoms */
+ rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i-1]+j]], dist);
+ length = norm(dist); /* determine distance between two atoms */
+ check_length(length, a[index[i-1]+j], a[index[i+1]+j]);
+
+ svmul(1.0/length, dist, Sz);
+ /* Sz is now the molecular axis Sz, normalized and all that */
+ }
+
+ /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
+ we can use the outer product of Cn-1->Cn and Cn+1->Cn, I hope */
+ rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], tmp1);
+ rvec_sub(x1[a[index[i-1]+j]], x1[a[index[i]+j]], tmp2);
+ cprod(tmp1, tmp2, Sx);
+ svmul(1.0/norm(Sx), Sx, Sx);
+
+ /* now we can get Sy from the outer product of Sx and Sz */
+ cprod(Sz, Sx, Sy);
+ svmul(1.0/norm(Sy), Sy, Sy);
+
+ /* the square of cosine of the angle between dist and the axis.
+ Using the innerproduct, but two of the three elements are zero
+ Determine the sum of the orderparameter of all atoms in group
+ */
+ if (use_unitvector)
+ {
+ cossum[XX] = gmx::square(iprod(Sx, direction)); /* this is allowed, since Sa is normalized */
+ cossum[YY] = gmx::square(iprod(Sy, direction));
+ cossum[ZZ] = gmx::square(iprod(Sz, direction));
+ }
+ else
+ {
+ cossum[XX] = gmx::square(Sx[axis]); /* this is allowed, since Sa is normalized */
+ cossum[YY] = gmx::square(Sy[axis]);
+ cossum[ZZ] = gmx::square(Sz[axis]);
+ }
+
+ for (m = 0; m < DIM; m++)
+ {
+ frameorder[m] += 0.5 * (3.0 * cossum[m] - 1.0);
+ }
+
+ if (bSliced)
+ {
+ /* get average coordinate in box length for slicing,
+ determine which slice atom is in, increase count for that
+ slice. slFrameorder and slOrder are reals, not
+ rvecs. Only the component [axis] of the order tensor is
+ kept, until I find it necessary to know the others too
+ */
+
+ z1 = x1[a[index[i-1]+j]][axis];
+ z2 = x1[a[index[i+1]+j]][axis];
+ z_ave = 0.5 * (z1 + z2);
- z_ave += box[axis][axis];
- }
- if (z_ave > box[axis][axis])
- {
- z_ave -= box[axis][axis];
++ slice = (int)((nslices*z_ave)/box[axis][axis]);
++ while (slice < 0)
+ {
- slice = static_cast<int>((0.5 + (z_ave / (*slWidth))) - 1);
++ slice += nslices;
+ }
++ slice = slice % nslices;
+
+ slCount[slice]++; /* determine slice, increase count */
+
+ slFrameorder[slice] += 0.5 * (3 * cossum[axis] - 1);
+ }
+ else if (permolecule)
+ {
+ /* store per-molecule order parameter
+ * To just track single-axis order: (*slOrder)[j][i] += 0.5 * (3 * iprod(cossum,direction) - 1);
+ * following is for Scd order: */
+ (*slOrder)[j][i] += -1* (1.0/3.0 * (3 * cossum[XX] - 1) + 1.0/3.0 * 0.5 * (3.0 * cossum[YY] - 1));
+ }
+ if (distcalc)
+ {
+ if (radial)
+ {
+ /* bin order parameter by arc distance from reference group*/
+ arcdist = gmx_angle(dvec, direction);
+ (*distvals)[j][i] += arcdist;
+ }
+ else if (i == 1)
+ {
+ /* Want minimum lateral distance to first group calculated */
+ tmpdist = trace(box); /* should be max value */
+ for (k = 0; k < distsize; k++)
+ {
+ pbc_dx(&pbc, x1[distidx[k]], x1[a[index[i]+j]], dvec);
+ /* at the moment, just remove dvec[axis] */
+ dvec[axis] = 0;
+ tmpdist = std::min(tmpdist, norm2(dvec));
+ }
+ //fprintf(stderr, "Min dist %f; trace %f\n", tmpdist, trace(box));
+ (*distvals)[j][i] += std::sqrt(tmpdist);
+ }
+ }
+ } /* end loop j, over all atoms in group */
+
+ for (m = 0; m < DIM; m++)
+ {
+ (*order)[i][m] += (frameorder[m]/size);
+ }
+
+ if (!permolecule)
+ { /*Skip following if doing per-molecule*/
+ for (k = 0; k < nslices; k++)
+ {
+ if (slCount[k]) /* if no elements, nothing has to be added */
+ {
+ (*slOrder)[k][i] += slFrameorder[k]/slCount[k];
+ slFrameorder[k] = 0; slCount[k] = 0;
+ }
+ }
+ } /* end loop i, over all groups in indexfile */
+ }
+ nr_frames++;
+
+ }
+ while (read_next_x(oenv, status, &t, x0, box));
+ /*********** done with status file **********/
+
+ fprintf(stderr, "\nRead trajectory. Printing parameters to file\n");
+ gmx_rmpbc_done(gpbc);
+
+ /* average over frames */
+ for (i = 1; i < ngrps - 1; i++)
+ {
+ svmul(1.0/nr_frames, (*order)[i], (*order)[i]);
+ fprintf(stderr, "Atom %d Tensor: x=%g , y=%g, z=%g\n", i, (*order)[i][XX],
+ (*order)[i][YY], (*order)[i][ZZ]);
+ if (bSliced || permolecule)
+ {
+ for (k = 0; k < nslices; k++)
+ {
+ (*slOrder)[k][i] /= nr_frames;
+ }
+ }
+ if (distcalc)
+ {
+ for (k = 0; k < nslices; k++)
+ {
+ (*distvals)[k][i] /= nr_frames;
+ }
+ }
+ }
+
+ if (bUnsat)
+ {
+ fprintf(stderr, "Average angle between double bond and normal: %f\n",
+ 180*sdbangle/(nr_frames * size*M_PI));
+ }
+
+ sfree(x0); /* free memory used by coordinate arrays */
+ sfree(x1);
+ if (comidx != NULL)
+ {
+ sfree(comidx);
+ }
+ if (distidx != NULL)
+ {
+ sfree(distidx);
+ }
+ if (grpname != NULL)
+ {
+ sfree(grpname);
+ }
+}
+
+
+void order_plot(rvec order[], real *slOrder[], const char *afile, const char *bfile,
+ const char *cfile, int ngrps, int nslices, real slWidth, gmx_bool bSzonly,
+ gmx_bool permolecule, real **distvals, const gmx_output_env_t *oenv)
+{
+ FILE *ord, *slOrd; /* xvgr files with order parameters */
+ int atom, slice; /* atom corresponding to order para.*/
+ char buf[256]; /* for xvgr title */
+ real S; /* order parameter averaged over all atoms */
+
+ if (permolecule)
+ {
+ sprintf(buf, "Scd order parameters");
+ ord = xvgropen(afile, buf, "Atom", "S", oenv);
+ sprintf(buf, "Orderparameters per atom per slice");
+ slOrd = xvgropen(bfile, buf, "Molecule", "S", oenv);
+ for (atom = 1; atom < ngrps - 1; atom++)
+ {
+ fprintf(ord, "%12d %12g\n", atom, -1.0 * (2.0/3.0 * order[atom][XX] +
+ 1.0/3.0 * order[atom][YY]));
+ }
+
+ for (slice = 0; slice < nslices; slice++)
+ {
+ fprintf(slOrd, "%12d\t", slice);
+ if (distvals)
+ {
+ fprintf(slOrd, "%12g\t", distvals[slice][1]); /*use distance value at second carbon*/
+ }
+ for (atom = 1; atom < ngrps - 1; atom++)
+ {
+ fprintf(slOrd, "%12g\t", slOrder[slice][atom]);
+ }
+ fprintf(slOrd, "\n");
+ }
+
+ }
+ else if (bSzonly)
+ {
+ sprintf(buf, "Orderparameters Sz per atom");
+ ord = xvgropen(afile, buf, "Atom", "S", oenv);
+ fprintf(stderr, "ngrps = %d, nslices = %d", ngrps, nslices);
+
+ sprintf(buf, "Orderparameters per atom per slice");
+ slOrd = xvgropen(bfile, buf, "Slice", "S", oenv);
+
+ for (atom = 1; atom < ngrps - 1; atom++)
+ {
+ fprintf(ord, "%12d %12g\n", atom, order[atom][ZZ]);
+ }
+
+ for (slice = 0; slice < nslices; slice++)
+ {
+ S = 0;
+ for (atom = 1; atom < ngrps - 1; atom++)
+ {
+ S += slOrder[slice][atom];
+ }
+ fprintf(slOrd, "%12g %12g\n", slice*slWidth, S/atom);
+ }
+
+ }
+ else
+ {
+ sprintf(buf, "Order tensor diagonal elements");
+ ord = xvgropen(afile, buf, "Atom", "S", oenv);
+ sprintf(buf, "Deuterium order parameters");
+ slOrd = xvgropen(cfile, buf, "Atom", "Scd", oenv);
+
+ for (atom = 1; atom < ngrps - 1; atom++)
+ {
+ fprintf(ord, "%12d %12g %12g %12g\n", atom, order[atom][XX],
+ order[atom][YY], order[atom][ZZ]);
+ fprintf(slOrd, "%12d %12g\n", atom, -1.0 * (2.0/3.0 * order[atom][XX] +
+ 1.0/3.0 * order[atom][YY]));
+ }
+
+ }
+ xvgrclose(ord);
+ xvgrclose(slOrd);
+}
+
+void write_bfactors(t_filenm *fnm, int nfile, int *index, int *a, int nslices, int ngrps, real **order, const t_topology *top, real **distvals, gmx_output_env_t *oenv)
+{
+ /*function to write order parameters as B factors in PDB file using
+ first frame of trajectory*/
+ t_trxstatus *status;
+ t_trxframe fr, frout;
+ t_atoms useatoms;
+ int i, j, ctr, nout;
+
+ ngrps -= 2; /*we don't have an order parameter for the first or
+ last atom in each chain*/
+ nout = nslices*ngrps;
+ read_first_frame(oenv, &status, ftp2fn(efTRX, nfile, fnm), &fr, TRX_NEED_X);
+
+ close_trj(status);
+ frout = fr;
+ frout.natoms = nout;
+ frout.bF = FALSE;
+ frout.bV = FALSE;
+ frout.x = 0;
+ snew(frout.x, nout);
+
+ init_t_atoms(&useatoms, nout, TRUE);
+ useatoms.nr = nout;
+
+ /*initialize PDBinfo*/
+ for (i = 0; i < useatoms.nr; ++i)
+ {
+ useatoms.pdbinfo[i].type = 0;
+ useatoms.pdbinfo[i].occup = 0.0;
+ useatoms.pdbinfo[i].bfac = 0.0;
+ useatoms.pdbinfo[i].bAnisotropic = FALSE;
+ }
+
+ for (j = 0, ctr = 0; j < nslices; j++)
+ {
+ for (i = 0; i < ngrps; i++, ctr++)
+ {
+ /*iterate along each chain*/
+ useatoms.pdbinfo[ctr].bfac = order[j][i+1];
+ if (distvals)
+ {
+ useatoms.pdbinfo[ctr].occup = distvals[j][i+1];
+ }
+ copy_rvec(fr.x[a[index[i+1]+j]], frout.x[ctr]);
+ useatoms.atomname[ctr] = top->atoms.atomname[a[index[i+1]+j]];
+ useatoms.atom[ctr] = top->atoms.atom[a[index[i+1]+j]];
+ useatoms.nres = std::max(useatoms.nres, useatoms.atom[ctr].resind+1);
+ useatoms.resinfo[useatoms.atom[ctr].resind] = top->atoms.resinfo[useatoms.atom[ctr].resind]; /*copy resinfo*/
+ }
+ }
+
+ write_sto_conf(opt2fn("-ob", nfile, fnm), "Order parameters", &useatoms, frout.x, NULL, frout.ePBC, frout.box);
+
+ sfree(frout.x);
+ done_atom(&useatoms);
+}
+
+int gmx_order(int argc, char *argv[])
+{
+ const char *desc[] = {
+ "[THISMODULE] computes the order parameter per atom for carbon tails. For atom i the",
+ "vector i-1, i+1 is used together with an axis. ",
+ "The index file should contain only the groups to be used for calculations,",
+ "with each group of equivalent carbons along the relevant acyl chain in its own",
+ "group. There should not be any generic groups (like System, Protein) in the index",
+ "file to avoid confusing the program (this is not relevant to tetrahedral order",
+ "parameters however, which only work for water anyway).[PAR]",
+ "[THISMODULE] can also give all",
+ "diagonal elements of the order tensor and even calculate the deuterium",
+ "order parameter Scd (default). If the option [TT]-szonly[tt] is given, only one",
+ "order tensor component (specified by the [TT]-d[tt] option) is given and the",
+ "order parameter per slice is calculated as well. If [TT]-szonly[tt] is not",
+ "selected, all diagonal elements and the deuterium order parameter is",
+ "given.[PAR]"
+ "The tetrahedrality order parameters can be determined",
+ "around an atom. Both angle an distance order parameters are calculated. See",
+ "P.-L. Chau and A.J. Hardwick, Mol. Phys., 93, (1998), 511-518.",
+ "for more details."
+ };
+
+ static int nslices = 1; /* nr of slices defined */
+ static gmx_bool bSzonly = FALSE; /* True if only Sz is wanted */
+ static gmx_bool bUnsat = FALSE; /* True if carbons are unsat. */
+ static const char *normal_axis[] = { NULL, "z", "x", "y", NULL };
+ static gmx_bool permolecule = FALSE; /*compute on a per-molecule basis */
+ static gmx_bool radial = FALSE; /*compute a radial membrane normal */
+ static gmx_bool distcalc = FALSE; /*calculate distance from a reference group */
+ t_pargs pa[] = {
+ { "-d", FALSE, etENUM, {normal_axis},
+ "Direction of the normal on the membrane" },
+ { "-sl", FALSE, etINT, {&nslices},
+ "Calculate order parameter as function of box length, dividing the box"
+ " into this number of slices." },
+ { "-szonly", FALSE, etBOOL, {&bSzonly},
+ "Only give Sz element of order tensor. (axis can be specified with [TT]-d[tt])" },
+ { "-unsat", FALSE, etBOOL, {&bUnsat},
+ "Calculate order parameters for unsaturated carbons. Note that this can"
+ "not be mixed with normal order parameters." },
+ { "-permolecule", FALSE, etBOOL, {&permolecule},
+ "Compute per-molecule Scd order parameters" },
+ { "-radial", FALSE, etBOOL, {&radial},
+ "Compute a radial membrane normal" },
+ { "-calcdist", FALSE, etBOOL, {&distcalc},
+ "Compute distance from a reference" },
+ };
+
+ rvec *order; /* order par. for each atom */
+ real **slOrder; /* same, per slice */
+ real slWidth = 0.0; /* width of a slice */
+ char **grpname; /* groupnames */
+ int ngrps, /* nr. of groups */
+ i,
+ axis = 0; /* normal axis */
+ t_topology *top; /* topology */
+ int ePBC;
+ int *index, /* indices for a */
+ *a; /* atom numbers in each group */
+ t_blocka *block; /* data from index file */
+ t_filenm fnm[] = { /* files for g_order */
+ { efTRX, "-f", NULL, ffREAD }, /* trajectory file */
+ { efNDX, "-n", NULL, ffREAD }, /* index file */
+ { efNDX, "-nr", NULL, ffOPTRD }, /* index for radial axis calculation */
+ { efTPR, NULL, NULL, ffREAD }, /* topology file */
+ { efXVG, "-o", "order", ffWRITE }, /* xvgr output file */
+ { efXVG, "-od", "deuter", ffWRITE }, /* xvgr output file */
+ { efPDB, "-ob", NULL, ffOPTWR }, /* write Scd as B factors to PDB if permolecule */
+ { efXVG, "-os", "sliced", ffWRITE }, /* xvgr output file */
+ { efXVG, "-Sg", "sg-ang", ffOPTWR }, /* xvgr output file */
+ { efXVG, "-Sk", "sk-dist", ffOPTWR }, /* xvgr output file */
+ { efXVG, "-Sgsl", "sg-ang-slice", ffOPTWR }, /* xvgr output file */
+ { efXVG, "-Sksl", "sk-dist-slice", ffOPTWR }, /* xvgr output file */
+ };
+ gmx_bool bSliced = FALSE; /* True if box is sliced */
+#define NFILE asize(fnm)
+ real **distvals = NULL;
+ const char *sgfnm, *skfnm, *ndxfnm, *tpsfnm, *trxfnm;
+ gmx_output_env_t *oenv;
+
+ if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME,
+ NFILE, fnm, asize(pa), pa, asize(desc), desc, 0, NULL, &oenv))
+ {
+ return 0;
+ }
+ if (nslices < 1)
+ {
+ gmx_fatal(FARGS, "Can not have nslices < 1");
+ }
+ sgfnm = opt2fn_null("-Sg", NFILE, fnm);
+ skfnm = opt2fn_null("-Sk", NFILE, fnm);
+ ndxfnm = opt2fn_null("-n", NFILE, fnm);
+ tpsfnm = ftp2fn(efTPR, NFILE, fnm);
+ trxfnm = ftp2fn(efTRX, NFILE, fnm);
+
+ /* Calculate axis */
+ GMX_RELEASE_ASSERT(normal_axis[0] != NULL, "Options inconsistency; normal_axis[0] is NULL");
+ if (std::strcmp(normal_axis[0], "x") == 0)
+ {
+ axis = XX;
+ }
+ else if (std::strcmp(normal_axis[0], "y") == 0)
+ {
+ axis = YY;
+ }
+ else if (std::strcmp(normal_axis[0], "z") == 0)
+ {
+ axis = ZZ;
+ }
+ else
+ {
+ gmx_fatal(FARGS, "Invalid axis, use x, y or z");
+ }
+
+ switch (axis)
+ {
+ case 0:
+ fprintf(stderr, "Taking x axis as normal to the membrane\n");
+ break;
+ case 1:
+ fprintf(stderr, "Taking y axis as normal to the membrane\n");
+ break;
+ case 2:
+ fprintf(stderr, "Taking z axis as normal to the membrane\n");
+ break;
+ }
+
+ /* tetraheder order parameter */
+ if (skfnm || sgfnm)
+ {
+ /* If either of theoptions is set we compute both */
+ sgfnm = opt2fn("-Sg", NFILE, fnm);
+ skfnm = opt2fn("-Sk", NFILE, fnm);
+ calc_tetra_order_parm(ndxfnm, tpsfnm, trxfnm, sgfnm, skfnm, nslices, axis,
+ opt2fn("-Sgsl", NFILE, fnm), opt2fn("-Sksl", NFILE, fnm),
+ oenv);
+ /* view xvgr files */
+ do_view(oenv, opt2fn("-Sg", NFILE, fnm), NULL);
+ do_view(oenv, opt2fn("-Sk", NFILE, fnm), NULL);
+ if (nslices > 1)
+ {
+ do_view(oenv, opt2fn("-Sgsl", NFILE, fnm), NULL);
+ do_view(oenv, opt2fn("-Sksl", NFILE, fnm), NULL);
+ }
+ }
+ else
+ {
+ /* tail order parameter */
+
+ if (nslices > 1)
+ {
+ bSliced = TRUE;
+ fprintf(stderr, "Dividing box in %d slices.\n\n", nslices);
+ }
+
+ if (bSzonly)
+ {
+ fprintf(stderr, "Only calculating Sz\n");
+ }
+ if (bUnsat)
+ {
+ fprintf(stderr, "Taking carbons as unsaturated!\n");
+ }
+
+ top = read_top(ftp2fn(efTPR, NFILE, fnm), &ePBC); /* read topology file */
+
+ block = init_index(ftp2fn(efNDX, NFILE, fnm), &grpname);
+ index = block->index; /* get indices from t_block block */
+ a = block->a; /* see block.h */
+ ngrps = block->nr;
+
+ if (permolecule)
+ {
+ nslices = index[1] - index[0]; /*I think this assumes contiguous lipids in topology*/
+ fprintf(stderr, "Calculating Scd order parameters for each of %d molecules\n", nslices);
+ }
+
+ if (radial)
+ {
+ fprintf(stderr, "Calculating radial distances\n");
+ if (!permolecule)
+ {
+ gmx_fatal(FARGS, "Cannot yet output radial distances without permolecule\n");
+ }
+ }
+
+ /* show atomtypes, to check if index file is correct */
+ print_types(index, a, ngrps, grpname, top);
+
+ calc_order(ftp2fn(efTRX, NFILE, fnm), index, a, &order,
+ &slOrder, &slWidth, nslices, bSliced, bUnsat,
+ top, ePBC, ngrps, axis, permolecule, radial, distcalc, opt2fn_null("-nr", NFILE, fnm), &distvals, oenv);
+
+ if (radial)
+ {
+ ngrps--; /*don't print the last group--was used for
+ center-of-mass determination*/
+
+ }
+ order_plot(order, slOrder, opt2fn("-o", NFILE, fnm), opt2fn("-os", NFILE, fnm),
+ opt2fn("-od", NFILE, fnm), ngrps, nslices, slWidth, bSzonly, permolecule, distvals, oenv);
+
+ if (opt2bSet("-ob", NFILE, fnm))
+ {
+ if (!permolecule)
+ {
+ fprintf(stderr,
+ "Won't write B-factors with averaged order parameters; use -permolecule\n");
+ }
+ else
+ {
+ write_bfactors(fnm, NFILE, index, a, nslices, ngrps, slOrder, top, distvals, oenv);
+ }
+ }
+
+
+ do_view(oenv, opt2fn("-o", NFILE, fnm), NULL); /* view xvgr file */
+ do_view(oenv, opt2fn("-os", NFILE, fnm), NULL); /* view xvgr file */
+ do_view(oenv, opt2fn("-od", NFILE, fnm), NULL); /* view xvgr file */
+ }
+
+ if (distvals != NULL)
+ {
+ for (i = 0; i < nslices; ++i)
+ {
+ sfree(distvals[i]);
+ }
+ sfree(distvals);
+ }
+
+ return 0;
+}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff