Convert manual from LaTeX to Markup
[gromacs.git] / docs / reference-manual / topologies.rst
blob6058d1f328d038ad7b009fdc0d40e0a7b7c56dc0
1 Topologies
2 ==========
4 Introduction
5 ------------
7 |Gromacs| must know on which atoms and combinations of atoms the various
8 contributions to the potential functions (see chapter :ref:`ff`) must act.
9 It must also know what parameters must be applied to the various
10 functions. All this is described in the *topology* file :ref:`top`, which
11 lists the *constant attributes* of each atom. There are many more atom
12 types than elements, but only atom types present in biological systems
13 are parameterized in the force field, plus some metals, ions and
14 silicon. The bonded and special interactions are determined by fixed
15 lists that are included in the topology file. Certain non-bonded
16 interactions must be excluded (first and second neighbors), as these are
17 already treated in bonded interactions. In addition, there are *dynamic
18 attributes* of atoms - their positions, velocities and forces. These do
19 not strictly belong to the molecular topology, and are stored in the
20 coordinate file :ref:`gro` (positions and velocities), or
21 trajectory file :ref:`trr` (positions, velocities, forces).
23 This chapter describes the setup of the topology file, the :ref:`top` file and
24 the database files: what the parameters stand for and how/where to
25 change them if needed. First, all file formats are explained. Section
26 :ref:`fffiles` describes the organization of the files in each force
27 field.
29 **Note:** if you construct your own topologies, we encourage you to
30 upload them to our topology archive at our `webpage`_! Just imagine how thankful
31 you’d have been if your topology had been available there before you
32 started. The same goes for new force fields or modified versions of the
33 standard force fields - contribute them to the force field archive!
35 .. _homepage: `webpage`_
37 Particle type
38 -------------
40 In |Gromacs|, there are three types of
41 particles
42 , see :numref:`Table %s <tab-ptype>`. Only regular atoms and virtual
43 interaction sites are used in |Gromacs|; shells are necessary for
44 polarizable models like the Shell-Water models \ :ref:`45 <refMaaren2001a>`.
46 .. _tab-ptype:
48 .. table:: Particle types in |Gromacs|
50            +--------------+----------+
51            | Particle     | Symbol   |
52            +==============+==========+
53            | atom         | A        |
54            +--------------+----------+
55            | shell        | S        |
56            +--------------+----------+
57            | virtual side | V (or D) |
58            +--------------+----------+
61 .. _atomtype:
63 Atom types
64 ~~~~~~~~~~
66 Each force field defines a set of atom
67 types,
68 which have a characteristic name or number, and mass (in a.m.u.). These
69 listings are found in the ``atomtypes.atp`` file (:ref:`atp` =
70 **a**\ tom **t**\ ype **p**\ arameter file). Therefore, it is in this
71 file that you can begin to change and/or add an atom type. A sample from
72 the ``gromos43a1.ff`` force field is listed below.
76      |  O  15.99940 ;     carbonyl oxygen (C=O)
77      | OM  15.99940 ;     carboxyl oxygen (CO-)
78      | OA  15.99940 ;     hydroxyl, sugar or ester oxygen
79      | OW  15.99940 ;     water oxygen
80      |  N  14.00670 ;     peptide nitrogen (N or NH)
81      | NT  14.00670 ;     terminal nitrogen (NH2)
82      | NL  14.00670 ;     terminal nitrogen (NH3)
83      | NR  14.00670 ;     aromatic nitrogen
84      | NZ  14.00670 ;     Arg NH (NH2)
85      | NE  14.00670 ;     Arg NE (NH)
86      |  C  12.01100 ;     bare carbon
87      |CH1  13.01900 ;     aliphatic or sugar CH-group
88      |CH2  14.02700 ;     aliphatic or sugar CH2-group
89      |CH3  15.03500 ;     aliphatic CH3-group
91 **Note:** |Gromacs| makes use of the atom types as a name, *not* as a
92 number (as *e.g.* in GROMOS).
94 .. _vsitetop:
96 Virtual sites
97 ~~~~~~~~~~~~~
99 Some force fields use virtual interaction sites (interaction sites that
100 are constructed from other particle positions) on which certain
101 interactions are located (*e.g.* on benzene rings, to reproduce the
102 correct quadrupole). This is described in sec. :ref:`virtualsites`.
104 To make virtual sites in your system, you should include a section
105 ``[ virtual_sites? ]`` (for backward compatibility the old
106 name ``[ dummies? ]`` can also be used) in your topology
107 file, where the ``?`` stands for the number constructing
108 particles for the virtual site. This will be `:ref:`2`` for
109 type 2, `:ref:`3`` for types 3, 3fd, 3fad and 3out and
110 `:ref:`4`` for type 4fdn. The last of these replace an older
111 4fd type (with the ‘type’ value 1) that could occasionally be unstable;
112 while it is still supported internally in the code, the old 4fd type
113 should not be used in new input files. The different types are explained
114 in sec. :ref:`virtualsites`.
116 Parameters for type 2 should look like this:
120     [ virtual_sites2 ]
121     ; Site  from        funct  a
122     5       1     2     1      0.7439756
124 for type 3 like this:
128     [ virtual_sites3 ]
129     ; Site  from               funct   a          b
130     5       1     2     3      1       0.7439756  0.128012
132 for type 3fd like this:
136     [ virtual_sites3 ]
137     ; Site  from               funct   a          d
138     5       1     2     3      2       0.5        -0.105
140 for type 3fad like this:
144     [ virtual_sites3 ]
145     ; Site  from               funct   theta      d
146     5       1     2     3      3       120        0.5
148 for type 3out like this:
152     [ virtual_sites3 ]
153     ; Site  from               funct   a          b          c
154     5       1     2     3      4       -0.4       -0.4       6.9281
156 for type 4fdn like this:
160     [ virtual_sites4 ]
161     ; Site  from                      funct   a          b          c
162     5       1     2     3     4       2       1.0        0.9       0.105
164 This will result in the construction of a virtual site, number 5 (first
165 column ``Site``), based on the positions of the atoms
166 whose indices are 1 and 2 or 1, 2 and 3 or 1, 2, 3 and 4 (next two,
167 three or four columns ``from``) following the rules
168 determined by the function number (next column ``funct``)
169 with the parameters specified (last one, two or three columns
170 ``a b . .``). Obviously, the atom numbers (including
171 virtual site number) depend on the molecule. It may be instructive to
172 study the topologies for TIP4P or TIP5P water models that are included
173 with the |Gromacs| distribution.
175 **Note** that if any constant bonded interactions are defined between
176 virtual sites and/or normal atoms, they will be removed by
177 :ref:`grompp <gmx grompp>` (unless the option ``-normvsbds`` is used). This
178 removal of bonded interactions is done after generating exclusions, as
179 the generation of exclusions is based on “chemically” bonded
180 interactions.
182 Virtual sites can be constructed in a more generic way using basic
183 geometric parameters. The directive that can be used is ``[ virtual_sitesn ]``. Required
184 parameters are listed in :numref:`Table %s <tab-topfile2>`. An example entry for
185 defining a virtual site at the center of geometry of a given set of
186 atoms might be:
190     [ virtual_sitesn ]
191     ; Site   funct    from
192     5        1        1     2     3     4
194 Parameter files
195 ---------------
197 Atoms
198 ~~~~~
200 The *static* properties (see  :numref:`Table %s <tab-statprop>`) assigned to the atom
201 types are assigned based on data in several places. The mass is listed
202 in ``atomtypes.atp`` (see :ref:`atomtype`), whereas the charge is listed
203 in :ref:`rtp` (:ref:`rtp` = **r**\ esidue **t**\ opology **p**\ arameter file,
204 see :ref:`rtp`). This implies that the charges are only defined in the
205 building blocks of amino acids, nucleic acids or otherwise, as defined
206 by the user. When generating a :ref:`topology <top>` using the :ref:`pdb2gmx <gmx pdb2gmx>`
207 program, the information from these files is combined.
209 .. _tab-statprop:
211 .. table:: Static atom type properties in |Gromacs|
213            +----------+------------------+----------+
214            | Property | Symbol           | Unit     |
215            +==========+==================+==========+
216            | Type     | -                | -        |
217            +----------+------------------+----------+
218            | Mass     | m                | a.m.u.   |
219            +----------+------------------+----------+
220            | Charge   | q                | electron |
221            +----------+------------------+----------+
222            | epsilon  | :math:`\epsilon` | kJ/mol   |
223            +----------+------------------+----------+
224            | sigma    | :math:`\sigma`   | nm       |
225            +----------+------------------+----------+
228 .. _nbpar:
230 Non-bonded parameters
231 ~~~~~~~~~~~~~~~~~~~~~
233 The non-bonded parameters consist of the van der Waals parameters V (``c6``
234 or :math:`\sigma`, depending on the combination rule) and W (``c12`` or
235 :math:`\epsilon`), as listed in the file ``ffnonbonded.itp``, where ``ptype`` is
236 the particle type (see :numref:`Table %s <tab-ptype>`). As with the bonded
237 parameters, entries in ``[ *type ]`` directives are applied to their counterparts in
238 the topology file. Missing parameters generate warnings, except as noted
239 below in section :ref:`pairinteractions`.
243     [ atomtypes ]
244     ;name   at.num      mass      charge   ptype         V(c6)        W(c12)
245         O        8  15.99940       0.000       A   0.22617E-02   0.74158E-06
246        OM        8  15.99940       0.000       A   0.22617E-02   0.74158E-06
247        .....
249     [ nonbond_params ]
250       ; i    j func       V(c6)        W(c12)
251         O    O    1 0.22617E-02   0.74158E-06
252         O   OA    1 0.22617E-02   0.13807E-05
253         .....
255 **Note** that most of the included force fields also include the ``at.num.``
256 column, but this same information is implied in the OPLS-AA ``bond_type``
257 column. The interpretation of the parameters V and W depends on the
258 combination rule that was chosen in the ``[ defaults ]`` section of the topology file
259 (see :ref:`topfile`):
261 .. math::
263    \begin{aligned}
264    \mbox{for combination rule 1}: & &
265    \begin{array}{llllll}
266      \mbox{V}_{ii} & = & C^{(6)}_{i}  & = & 4\,\epsilon_i\sigma_i^{6} &
267      \mbox{[ kJ mol$^{-1}$ nm$^{6}$ ]}\\
268      \mbox{W}_{ii} & = & C^{(12)}_{i} & = & 4\,\epsilon_i\sigma_i^{12} &
269      \mbox{[ kJ mol$^{-1}$ nm$^{12}$ ]}\\
270    \end{array}
271    \\
272    \mbox{for combination rules 2 and 3}: & &
273    \begin{array}{llll}
274      \mbox{V}_{ii} & = & \sigma_i   & \mbox{[ nm ]} \\
275      \mbox{W}_{ii} & = & \epsilon_i & \mbox{[ kJ mol$^{-1}$ ]}
276    \end{array}\end{aligned}
278 Some or all combinations for different atom types can be given in the
279 ``[ nonbond_params ]`` section, again with parameters V and
280 W as defined above. Any combination that is not given will be computed
281 from the parameters for the corresponding atom types, according to the
282 combination rule:
284 .. math::
286    \begin{aligned}
287    \mbox{for combination rules 1 and 3}: & &
288    \begin{array}{lll}
289      C^{(6)}_{ij}  & = & \left(C^{(6)}_i\,C^{(6)}_j\right)^{\frac{1}{2}} \\
290      C^{(12)}_{ij} & = & \left(C^{(12)}_i\,C^{(12)}_j\right)^{\frac{1}{2}}
291    \end{array}
292    \\
293    \mbox{for combination rule 2}: & &
294    \begin{array}{lll}
295      \sigma_{ij}   & = & \frac{1}{2}(\sigma_i+\sigma_j) \\
296      \epsilon_{ij} & = & \sqrt{\epsilon_i\,\epsilon_j}
297    \end{array}\end{aligned}
299 When :math:`\sigma` and :math:`\epsilon` need to be supplied (rules 2
300 and 3), it would seem it is impossible to have a non-zero :math:`C^{12}`
301 combined with a zero :math:`C^6` parameter. However, providing a
302 negative :math:`\sigma` will do exactly that, such that :math:`C^6` is
303 set to zero and :math:`C^{12}` is calculated normally. This situation
304 represents a special case in reading the value of :math:`\sigma`, and
305 nothing more.
307 There is only one set of combination rules for Buckingham potentials:
309 .. math::
311    \begin{array}{rcl}
312    A_{ij}   &=& \left(A_{ii} \, A_{jj}\right)^{1/2}    \\
313    B_{ij}   &=& 2 / \left(\frac{1}{B_{ii}} + \frac{1}{B_{jj}}\right)        \\
314    C_{ij}   &=& \left(C_{ii} \, C_{jj}\right)^{1/2}
315    \end{array}
317 Bonded parameters
318 ~~~~~~~~~~~~~~~~~
320 The bonded
321 parameters
322 (*i.e.* bonds, bond angles, improper and proper dihedrals) are listed in
323 ``ffbonded.itp``.  The entries in this database describe,
324 respectively, the atom types in the interactions, the type of the
325 interaction, and the parameters associated with that interaction. These
326 parameters are then read by
327 :ref:`grompp <gmx grompp>` when processing a
328 topology and applied to the relevant bonded parameters, *i.e.*
329 ``bondtypes`` are applied to entries in the
330 ``[ bonds ]`` directive, etc. Any bonded parameter that is
331 missing from the relevant :``[ *type ]`` directive generates
332 a fatal error. The types of interactions are listed in
333 :numref:`Table %s <tab-topfile2>`. Example excerpts from such files
334 follow:
338     [ bondtypes ]
339       ; i    j func        b0          kb
340         C    O    1   0.12300     502080.
341         C   OM    1   0.12500     418400.
342         ......
344     [ angletypes ]
345       ; i    j    k func       th0         cth
346        HO   OA    C    1   109.500     397.480
347        HO   OA  CH1    1   109.500     397.480
348        ......
350     [ dihedraltypes ]
351       ; i    l func        q0          cq
352      NR5*  NR5    2     0.000     167.360
353      NR5* NR5*    2     0.000     167.360
354      ......
356     [ dihedraltypes ]
357       ; j    k func      phi0          cp   mult
358         C   OA    1   180.000      16.736      2
359         C    N    1   180.000      33.472      2
360         ......
362     [ dihedraltypes ]
363     ;
364     ; Ryckaert-Bellemans Dihedrals
365     ;
366     ; aj    ak      funct
367     CP2     CP2     3       9.2789  12.156  -13.120 -3.0597 26.240  -31.495
369 In the ``ffbonded.itp`` file, you can add bonded parameters.
370 If you want to include parameters for new atom types, make sure you
371 define them in ``atomtypes.atp`` as well.
373 For most interaction types, bonded parameters are searched and assigned
374 using an exact match for all type names and allowing only a single set
375 of parameters. The exception to this rule are
376 dihedral
377 parameters. For
378 ``[ dihedraltypes ]`` wildcard atom type names can be
379 specified with the letter ``X`` in one or more of the four
380 positions. Thus one can for example assign proper dihedral parameters
381 based on the types of the middle two atoms. The parameters for the entry
382 with the most exact matches, i.e. the least wildcard matches, will be
383 used. Note that |Gromacs| versions older than 5.1.3 used the first match,
384 which means that a full match would be ignored if it is preceded by an
385 entry that matches on wildcards. Thus it is suggested to put wildcard
386 entries at the end, in case someone might use a forcefield with older
387 versions of |Gromacs|. In addition there is a dihedral type 9 which adds
388 the possibility of assigning multiple dihedral potentials, useful for
389 combining terms with different multiplicities. The different dihedral
390 potential parameter sets should be on directly adjacent lines in the
391 ``[ dihedraltypes ]`` section.
393 Molecule definition
394 -------------------
396 Moleculetype entries
397 ~~~~~~~~~~~~~~~~~~~~
399 An organizational structure that usually corresponds to molecules is the
400 ``[ moleculetype ]`` entry. This entry serves two main
401 purposes. One is to give structure to the topology file(s), usually
402 corresponding to real molecules. This makes the topology easier to read
403 and writing it less labor intensive. A second purpose is computational
404 efficiency. The system definition that is kept in memory is proportional
405 in size of the ``moleculetype`` definitions. If a molecule
406 is present in 100000 copies, this saves a factor of 100000 in memory,
407 which means the system usually fits in cache, which can improve
408 performance tremendously. Interactions that correspond to chemical
409 bonds, that generate exclusions, can only be defined between atoms
410 within a ``moleculetype``. It is allowed to have multiple
411 molecules which are not covalently bonded in one
412 ``moleculetype`` definition. Molecules can be made
413 infinitely long by connecting to themselves over periodic boundaries.
414 When such periodic molecules are present, an option in the
415 :ref:`mdp` file needs to be set to tell |Gromacs| not to attempt
416 to make molecules that are broken over periodic boundaries whole again.
418 Intermolecular interactions
419 ~~~~~~~~~~~~~~~~~~~~~~~~~~~
421 In some cases, one would like atoms in different molecules to also
422 interact with other interactions than the usual non-bonded interactions.
423 This is often the case in binding studies. When the molecules are
424 covalently bound, e.g. a ligand binding covalently to a protein, they
425 are effectively one molecule and they should be defined in one
426 ``[ moleculetype ]`` entry. Note that
427 :ref:`pdb2gmx <gmx pdb2gmx>` has an option to put two or more molecules in
428 one ``[ moleculetype ]`` entry. When molecules are not
429 covalently bound, it is much more convenient to use separate
430 ``moleculetype`` definitions and specify the intermolecular
431 interactions in the ``[ intermolecular_interactions]``
432 section. In this section, which is placed at the end of the topology
433 (see :numref:`Table %s <tab-topfile1>`), normal bonded interactions
434 can be specified using global atom indices. The only restrictions are
435 that no interactions can be used that generates exclusions and no
436 constraints can be used.
438 .. _pairinteractions:
440 Intramolecular pair interactions
441 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
443 Extra Lennard-Jones and electrostatic interactions between pairs of
444 atoms in a molecule can be added in the ``[ pairs ]`` section of a molecule
445 definition. The parameters for these interactions can be set
446 independently from the non-bonded interaction parameters. In the GROMOS
447 force fields, pairs are only used to modify the 1-4 interactions
448 (interactions of atoms separated by three bonds). In these force fields
449 the 1-4 interactions are excluded from the non-bonded interactions (see
450 sec. :ref:`excl`).
455     [ pairtypes ]
456       ; i    j func         cs6          cs12 ; THESE ARE 1-4 INTERACTIONS
457         O    O    1 0.22617E-02   0.74158E-06
458         O   OM    1 0.22617E-02   0.74158E-06
459         .....
461 The pair interaction parameters for the atom types in ``ffnonbonded.itp``
462 are listed in the ``[ pairtypes ]`` section. The GROMOS force fields list all these
463 interaction parameters explicitly, but this section might be empty for
464 force fields like OPLS that calculate the 1-4 interactions by uniformly
465 scaling the parameters. Pair parameters that are not present in the ``[ pairtypes ]``
466 section are only generated when ``gen-pairs`` is set to ``yes`` in the
467 ``[ defaults ]`` directive of ``forcefield.itp`` (see :ref:`topfile`). When ``gen-pairs`` is
468 set to ``no``, :ref:`grompp <gmx grompp>` will give a warning for each pair type for which no
469 parameters are given.
471 The normal pair interactions, intended for 1-4 interactions, have
472 function type 1. Function type 2 and the ``[ pairs_nb ]`` are intended for free-energy
473 simulations. When determining hydration free energies, the solute needs
474 to be decoupled from the solvent. This can be done by adding a B-state
475 topology (see sec. :ref:`fecalc`) that uses zero for all solute
476 non-bonded parameters, *i.e.* charges and LJ parameters. However, the
477 free energy difference between the A and B states is not the total
478 hydration free energy. One has to add the free energy for reintroducing
479 the internal Coulomb and LJ interactions in the solute when in vacuum.
480 This second step can be combined with the first step when the Coulomb
481 and LJ interactions within the solute are not modified. For this
482 purpose, there is a pairs function type 2, which is identical to
483 function type 1, except that the B-state parameters are always identical
484 to the A-state parameters. For searching the parameters in the ``[ pairtypes ]`` section,
485 no distinction is made between function type 1 and 2. The pairs section
486 ``[ pairs_nb ]`` is intended to replace the non-bonded interaction. It uses the unscaled
487 charges and the non-bonded LJ parameters; it also only uses the A-state
488 parameters. **Note** that one should add exclusions for all atom pairs
489 listed in ``[ pairs_nb ]``, otherwise such pairs will also end up in the normal neighbor
490 lists.
492 Alternatively, this same behavior can be achieved without ever touching
493 the topology, by using the ``couple-moltype``, ``couple-lambda0``,
494 ``couple-lambda1``, and ``couple-intramol`` keywords. See sections
495 sec. :ref:`fecalc` and sec. :ref:`dgimplement` for more information.
497 All three pair types always use plain Coulomb interactions, even when
498 Reaction-field, PME, Ewald or shifted Coulomb interactions are selected
499 for the non-bonded interactions. Energies for types 1 and 2 are written
500 to the energy and log file in separate “LJ-14” and “Coulomb-14” entries
501 per energy group pair. Energies for ``[ pairs_nb ]`` are added to the “LJ-(SR)” and
502 “Coulomb-(SR)” terms.
504 .. _excl:
506 Exclusions
507 ~~~~~~~~~~
509 The exclusions for non-bonded interactions are generated by :ref:`grompp <gmx grompp>` for
510 neighboring atoms up to a certain number of bonds away, as defined in
511 the ``[ moleculetype ]`` section in the topology file (see :ref:`topfile`). Particles are
512 considered bonded when they are connected by “chemical” bonds (``[ bonds ]`` types 1
513 to 5, 7 or 8) or constraints (``[ constraints ]`` type 1). Type 5 ``[ bonds ]`` can be used to create a
514 connection between two atoms without creating an interaction. There is a
515 harmonic interaction (``[ bonds ]`` type 6) that does not connect the atoms by a
516 chemical bond. There is also a second constraint type (``[ constraints ]`` type 2) that
517 fixes the distance, but does not connect the atoms by a chemical bond.
518 For a complete list of all these interactions, see :numref:`Table %s <tab-topfile2>`.
520 Extra exclusions within a molecule can be added manually in a
521 ``[ exclusions ]`` section. Each line should start with one
522 atom index, followed by one or more atom indices. All non-bonded
523 interactions between the first atom and the other atoms will be
524 excluded.
526 When all non-bonded interactions within or between groups of atoms need
527 to be excluded, is it more convenient and much more efficient to use
528 energy monitor group exclusions (see sec. :ref:`groupconcept`).
530 .. _constraintalg:
532 Constraint algorithms
533 ---------------------
535 Constraints are defined in the ``[ constraints ]`` section. The format is two atom numbers
536 followed by the function type, which can be 1 or 2, and the constraint
537 distance. The only difference between the two types is that type 1 is
538 used for generating exclusions and type 2 is not (see sec. :ref:`excl`).
539 The distances are constrained using the LINCS or the SHAKE algorithm,
540 which can be selected in the :ref:`mdp` file. Both types of constraints can be
541 perturbed in free-energy calculations by adding a second constraint
542 distance (see :ref:`constraintforce`). Several types of bonds and
543 angles (see :numref:`Table %s <tab-topfile2>`) can be converted automatically to
544 constraints by :ref:`grompp <gmx grompp>`. There are several options for this in the :ref:`mdp`
545 file.
547 We have also implemented the SETTLE
548 algorithm \ :ref:`47 <refMiyamoto92>`, which is an analytical solution of SHAKE, specifically for
549 water. SETTLE can be selected in the topology file. See, for instance,
550 the SPC molecule definition:
554     [ moleculetype ]
555     ; molname       nrexcl
556     SOL             1
558     [ atoms ]
559     ; nr    at type res nr  ren nm  at nm   cg nr   charge
560     1       OW      1       SOL     OW1     1       -0.82
561     2       HW      1       SOL     HW2     1        0.41
562     3       HW      1       SOL     HW3     1        0.41
564     [ settles ]
565     ; OW    funct   doh     dhh
566     1       1       0.1     0.16333
568     [ exclusions ]
569     1       2       3
570     2       1       3
571     3       1       2
573 The ``[ settles ]`` directive defines the first atom of the
574 water molecule. The settle funct is always 1, and the distance between
575 O-H and H-H distances must be given. **Note** that the algorithm can
576 also be used for TIP3P and TIP4P \ :ref:`128 <refJorgensen83>`. TIP3P just has
577 another geometry. TIP4P has a virtual site, but since that is generated
578 it does not need to be shaken (nor stirred).
580 .. _pdb2gmxfiles:
582 :ref:`pdb2gmx <gmx pdb2gmx>` input files
583 ----------------------------------------
585 The |Gromacs| program :ref:`pdb2gmx <gmx pdb2gmx>` generates a topology for the input
586 coordinate file. Several formats are supported for that coordinate file,
587 but :ref:`pdb` is the most commonly-used format (hence the name :ref:`pdb2gmx <gmx pdb2gmx>`).
588 :ref:`pdb2gmx <gmx pdb2gmx>` searches for force fields in sub-directories of the |Gromacs|
589 ``share/top`` directory and your working directory. Force fields are
590 recognized from the file ``forcefield.itp`` in a directory with the
591 extension ``.ff``. The file ``forcefield.doc`` may be present, and if so, its
592 first line will be used by :ref:`pdb2gmx <gmx pdb2gmx>` to present a short description to the
593 user to help in choosing a force field. Otherwise, the user can choose a
594 force field with the ``-ff xxx`` command-line argument to :ref:`pdb2gmx <gmx pdb2gmx>`, which
595 indicates that a force field in a ``xxx.ff`` directory is desired. :ref:`pdb2gmx <gmx pdb2gmx>`
596 will search first in the working directory, then in the |Gromacs|
597 ``share/top`` directory, and use the first matching ``xxx.ff`` directory found.
599 Two general files are read by :ref:`pdb2gmx <gmx pdb2gmx>`: an atom type file (extension
600 :ref:`atp`, see :ref:`atomtype`) from the force-field directory, and a file
601 called ``residuetypes.dat`` from either the working directory, or the
602 |Gromacs| ``share/top`` directory. ``residuetypes.dat`` determines which residue
603 names are considered protein, DNA, RNA, water, and ions.
605 :ref:`pdb2gmx <gmx pdb2gmx>` can read one or multiple databases with topological information
606 for different types of molecules. A set of files belonging to one
607 database should have the same basename, preferably telling something
608 about the type of molecules (*e.g.* aminoacids, rna, dna). The possible
609 files are:
611 -  ``<basename>.rtp``
613 -  ``<basename>.r2b (optional)``
615 -  ``<basename>.arn (optional)``
617 -  ``<basename>.hdb (optional)``
619 -  ``<basename>.n.tdb (optional)``
621 -  ``<basename>.c.tdb (optional)``
623 Only the :ref:`rtp` file, which contains the topologies of the building
624 blocks, is mandatory. Information from other files will only be used for
625 building blocks that come from an :ref:`rtp` file with the same base name. The
626 user can add building blocks to a force field by having additional files
627 with the same base name in their working directory. By default, only
628 extra building blocks can be defined, but calling :ref:`pdb2gmx <gmx pdb2gmx>` with the ``-rtpo``
629 option will allow building blocks in a local file to replace the default
630 ones in the force field.
632 Residue database
633 ~~~~~~~~~~~~~~~~
635 The files holding the residue databases have the extension :ref:`rtp`.
636 Originally this file contained building blocks (amino acids) for
637 proteins, and is the |Gromacs| interpretation of the ``rt37c4.dat`` file of
638 GROMOS. So the residue database file contains information (bonds,
639 charges, charge groups, and improper dihedrals) for a frequently-used
640 building block. It is better *not* to change this file because it is
641 standard input for :ref:`pdb2gmx <gmx pdb2gmx>`, but if changes are needed make them in the
642 :ref:`top` file (see :ref:`topfile`), or in a :ref:`rtp` file in the working
643 directory as explained in sec. :ref:`pdb2gmxfiles`. Defining topologies
644 of new small molecules is probably easier by writing an include topology
645 file :ref:`itp` directly. This will be discussed in section :ref:`molitp`.
646 When adding a new protein residue to the database, don’t forget to add
647 the residue name to the residuetypes.dat file, so that :ref:`grompp <gmx grompp>`, :ref:`make_ndx <gmx make_ndx>`
648 and analysis tools can recognize the residue as a protein residue (see
649 :ref:`defaultgroups`).
651 The :ref:`rtp` files are only used by :ref:`pdb2gmx <gmx pdb2gmx>`. As mentioned before, the only
652 extra information this program needs from the :ref:`rtp` database is bonds,
653 charges of atoms, charge groups, and improper dihedrals, because the
654 rest is read from the coordinate input file. Some proteins contain
655 residues that are not standard, but are listed in the coordinate file.
656 You have to construct a building block for this “strange” residue,
657 otherwise you will not obtain a :ref:`top` file. This also holds for molecules
658 in the coordinate file such as ligands, polyatomic ions, crystallization
659 co-solvents, etc. The residue database is constructed in the following
660 way:
664     [ bondedtypes ]  ; mandatory
665     ; bonds  angles  dihedrals  impropers
666          1       1          1          2  ; mandatory
668     [ GLY ]  ; mandatory
670      [ atoms ]  ; mandatory 
671     ; name  type  charge  chargegroup 
672          N     N  -0.280     0
673          H     H   0.280     0
674         CA   CH2   0.000     1
675          C     C   0.380     2
676          O     O  -0.380     2
678      [ bonds ]  ; optional
679     ;atom1 atom2      b0      kb
680          N     H
681          N    CA
682         CA     C
683          C     O
684         -C     N
686      [ exclusions ]  ; optional
687     ;atom1 atom2
689      [ angles ]  ; optional
690     ;atom1 atom2 atom3    th0    cth
692      [ dihedrals ]  ; optional
693     ;atom1 atom2 atom3 atom4   phi0     cp   mult
695      [ impropers ]  ; optional
696     ;atom1 atom2 atom3 atom4     q0     cq
697          N    -C    CA     H
698         -C   -CA     N    -O
700     [ ZN ]
702      [ atoms ]
703         ZN    ZN   2.000     0
705 The file is free format; the only restriction is that there can be at
706 most one entry on a line. The first field in the file is the ``[ bondedtypes ]`` field,
707 which is followed by four numbers, indicating the interaction type for
708 bonds, angles, dihedrals, and improper dihedrals. The file contains
709 residue entries, which consist of atoms and (optionally) bonds, angles,
710 dihedrals, and impropers. The charge group codes denote the charge group
711 numbers. Atoms in the same charge group should always be ordered
712 consecutively. When using the hydrogen database with :ref:`pdb2gmx <gmx pdb2gmx>` for adding
713 missing hydrogens (see :ref:`hdb`), the atom names defined in the :ref:`rtp`
714 entry should correspond exactly to the naming convention used in the
715 hydrogen database. The atom names in the bonded interaction can be
716 preceded by a minus or a plus, indicating that the atom is in the
717 preceding or following residue respectively. Explicit parameters added
718 to bonds, angles, dihedrals, and impropers override the standard
719 parameters in the :ref:`itp` files. This should only be used in special cases.
720 Instead of parameters, a string can be added for each bonded
721 interaction. This is used in GROMOS-96 :ref:`rtp` files. These strings are
722 copied to the topology file and can be replaced by force-field
723 parameters by the C-preprocessor in :ref:`grompp <gmx grompp>` using ``#define`` statements.
725 :ref:`pdb2gmx <gmx pdb2gmx>` automatically generates all angles. This means
726 that for most force fields the ``[ angles ]`` field is only
727 useful for overriding :ref:`itp` parameters. For the GROMOS-96
728 force field the interaction number of all angles needs to be specified.
730 :ref:`pdb2gmx <gmx pdb2gmx>` automatically generates one proper dihedral for every rotatable
731 bond, preferably on heavy atoms. When the ``[ dihedrals ]`` field is used, no other
732 dihedrals will be generated for the bonds corresponding to the specified
733 dihedrals. It is possible to put more than one dihedral function on a
734 rotatable bond. In the case of CHARMM27 FF :ref:`pdb2gmx <gmx pdb2gmx>` can add correction
735 maps to the dihedrals using the default ``-cmap`` option. Please refer to
736 :ref:`charmmff` for more information.
738 :ref:`pdb2gmx <gmx pdb2gmx>` sets the number of exclusions to 3, which means
739 that interactions between atoms connected by at most 3 bonds are
740 excluded. Pair interactions are generated for all pairs of atoms that
741 are separated by 3 bonds (except pairs of hydrogens). When more
742 interactions need to be excluded, or some pair interactions should not
743 be generated, an ``[ exclusions ]`` field can be added,
744 followed by pairs of atom names on separate lines. All non-bonded and
745 pair interactions between these atoms will be excluded.
747 Residue to building block database
748 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
750 Each force field has its own naming convention for residues. Most
751 residues have consistent naming, but some, especially those with
752 different protonation states, can have many different names. The
753 :ref:`r2b` files are used to convert standard residue names to
754 the force-field build block names. If no :ref:`r2b` is present
755 in the force-field directory or a residue is not listed, the building
756 block name is assumed to be identical to the residue name. The
757 :ref:`r2b` can contain 2 or 5 columns. The 2-column format has
758 the residue name in the first column and the building block name in the
759 second. The 5-column format has 3 additional columns with the building
760 block for the residue occurring in the N-terminus, C-terminus and both
761 termini at the same time (single residue molecule). This is useful for,
762 for instance, the AMBER force fields. If one or more of the terminal
763 versions are not present, a dash should be entered in the corresponding
764 column.
766 There is a |Gromacs| naming convention for residues which is only apparent
767 (except for the :ref:`pdb2gmx <gmx pdb2gmx>` code) through the
768 :ref:`r2b` file and ``specbond.dat`` files. This
769 convention is only of importance when you are adding residue types to an
770 :ref:`rtp` file. The convention is listed in :numref:`Table %s <tab-r2b>`.
771 For special bonds with, for instance,
772 a heme group, the |Gromacs| naming convention is introduced through
773 ``specbond.dat`` (see :ref:`specbond`),
774 which can subsequently be translated by the :ref:`r2b` file,
775 if required.
777 .. |NDEL| replace:: N\ :math:`_\delta`
778 .. |NEPS| replace:: N\ :math:`_\epsilon`
780 .. _tab-r2b:
782 .. table:: Internal |Gromacs| residue naming convention.
784            +--------------+-----------------------------------------------------------+
785            | |Gromacs| ID | Residue                                                   |
786            +==============+===========================================================+
787            | ARG          | protonated arginine                                       |
788            +--------------+-----------------------------------------------------------+
789            | ARGN         | neutral arginine                                          |
790            +--------------+-----------------------------------------------------------+
791            | ASP          | negatively charged aspartic acid                          |
792            +--------------+-----------------------------------------------------------+
793            | ASPH         | neutral aspartic acid                                     |
794            +--------------+-----------------------------------------------------------+
795            | CYS          | neutral cysteine                                          |
796            +--------------+-----------------------------------------------------------+
797            | CYS2         | cysteine with sulfur bound to another cysteine or a heme  |
798            +--------------+-----------------------------------------------------------+
799            | GLU          |  negatively charged glutamic acid                         |
800            +--------------+-----------------------------------------------------------+
801            | GLUH         |  neutral glutamic acid                                    |
802            +--------------+------------------------------+----------------------------+
803            | HISD         | neutral histidine with |NDEL| protonated                  |
804            +--------------+-----------------------------------------------------------+
805            | HISE         | neutral histidine with |NEPS| protonated                  |
806            +--------------+------------------------------+----------------------------+
807            | HISH         | positive histidine with both |NDEL| and |NEPS| protonated |
808            +--------------+-----------------------------------------------------------+
809            | HIS1         | histidine bound to a heme                                 |
810            +--------------+-----------------------------------------------------------+
811            | LYSN         | neutral lysine                                            |
812            +--------------+-----------------------------------------------------------+
813            | LYS          | protonated lysine                                         |
814            +--------------+-----------------------------------------------------------+
815            | HEME         | heme                                                      |
816            +--------------+-----------------------------------------------------------+
819 Atom renaming database
820 ~~~~~~~~~~~~~~~~~~~~~~
822 Force fields often use atom names that do not follow IUPAC or PDB
823 convention. The :ref:`arn` database is used to translate the
824 atom names in the coordinate file to the force-field names. Atoms that
825 are not listed keep their names. The file has three columns: the
826 building block name, the old atom name, and the new atom name,
827 respectively. The residue name supports question-mark wildcards that
828 match a single character.
830 An additional general atom renaming file called
831 ``xlateat.dat`` is present in the ``share/top``
832 directory, which translates common non-standard atom names in the
833 coordinate file to IUPAC/PDB convention. Thus, when writing force-field
834 files, you can assume standard atom names and no further atom name
835 translation is required, except for translating from standard atom names
836 to the force-field ones.
838 Hydrogen database
839 ~~~~~~~~~~~~~~~~~
841 The hydrogen database is stored in :ref:`hdb` files. It contains information
842 for the :ref:`pdb2gmx <gmx pdb2gmx>` program on how to connect hydrogen atoms to existing
843 atoms. In versions of the database before |Gromacs| 3.3, hydrogen atoms
844 were named after the atom they are connected to: the first letter of the
845 atom name was replaced by an ‘H.’ In the versions from 3.3 onwards, the
846 H atom has to be listed explicitly, because the old behavior was
847 protein-specific and hence could not be generalized to other molecules.
848 If more than one hydrogen atom is connected to the same atom, a number
849 will be added to the end of the hydrogen atom name. For example, adding
850 two hydrogen atoms to ``ND2`` (in asparagine), the hydrogen atoms will
851 be named ``HD21`` and ``HD22``. This is important since atom naming in
852 the :ref:`rtp` file (see :ref:`rtp`) must be the same. The format of the
853 hydrogen database is as follows:
857     ; res   # additions
858             # H add type    H       i       j       k
859     ALA     1
860             1       1       H       N       -C      CA
861     ARG     4
862             1       2       H       N       CA      C
863             1       1       HE      NE      CD      CZ
864             2       3       HH1     NH1     CZ      NE
865             2       3       HH2     NH2     CZ      NE
867 On the first line we see the residue name (ALA or ARG) and the number of
868 kinds of hydrogen atoms that may be added to this residue by the
869 hydrogen database. After that follows one line for each addition, on
870 which we see:
872 -  The number of H atoms added
874 -  The method for adding H atoms, which can be any of:
876    #. | *one planar hydrogen*, *e.g.* *rings or peptide bond*
877       | One hydrogen atom (n) is generated, lying in the plane of atoms
878         (i,j,k) on the plane bisecting angle (j-i-k) at a distance of
879         0.1 nm from atom i, such that the angles (n-i-j) and (n-i-k) are
880         :math:`>` 90\ :math:`^{\rm o}`.
882    #. | *one single hydrogen*, *e.g.* *hydroxyl*
883       | One hydrogen atom (n) is generated at a distance of 0.1 nm from
884         atom i, such that angle (n-i-j)=109.5 degrees and dihedral
885         (n-i-j-k)=trans.
887    #. | *two planar hydrogens*, *e.g.* *ethylene -C=CH*:math:`_2`, *or amide
888         -C(=O)NH*:math:`_2`
889       | Two hydrogens (n1,n2) are generated at a distance of 0.1 nm from
890         atom i, such that angle (n1-i-j)=(n2-i-j)=120 degrees and
891         dihedral (n1-i-j-k)=cis and (n2-i-j-k)=trans, such that names
892         are according to IUPAC standards \ :ref:`129 <refiupac70>`.
894    #. | *two or three tetrahedral hydrogens*, *e.g.* *-CH*:math:`_3`
895       | Three (n1,n2,n3) or two (n1,n2) hydrogens are generated at a
896         distance of 0.1 nm from atom i, such that angle
897         (n1-i-j)=(n2-i-j)=(n3-i-j)=109.47:math:`^{\rm o}`, dihedral
898         (n1-i-j-k)=trans, (n2-i-j-k)=trans+120 and
899         (n3-i-j-k)=trans+240:math:`^{\rm o}`.
901    #. | *one tetrahedral hydrogen*, *e.g.* *C*\ :math:`_3`\* CH*
902       | One hydrogen atom (n:math:`^\prime`) is generated at a distance
903         of 0.1 nm from atom i in tetrahedral conformation such that
904         angle
905         (n:math:`^\prime`-i-j)=(n:math:`^\prime`-i-k)=(n:math:`^\prime`-i-l)=109.47:math:`^{\rm o}`.
907    #. | *two tetrahedral hydrogens*, *e.g.* *C-CH*\ :math:`_2`\*-C*
908       | Two hydrogen atoms (n1,n2) are generated at a distance of 0.1 nm
909         from atom i in tetrahedral conformation on the plane bisecting
910         angle j-i-k with angle
911         (n1-i-n2)=(n1-i-j)=(n1-i-k)=109.47:math:`^{\rm o}`.
913    #. | *two water hydrogens*
914       | Two hydrogens are generated around atom i according to
915         SPC \ :ref:`80 <refBerendsen81>` water geometry. The symmetry
916         axis will alternate between three coordinate axes in both
917         directions.
919    #. | *three water “hydrogens”*
920       | Two hydrogens are generated around atom i according to
921         SPC \ :ref:`80 <refBerendsen81>` water geometry. The symmetry
922         axis will alternate between three coordinate axes in both
923         directions. In addition, an extra particle is generated on the
924         position of the oxygen with the first letter of the name
925         replaced by ‘M’. This is for use with four-atom water models
926         such as TIP4P \ :ref:`128 <refJorgensen83>`.
928    #. | *four water “hydrogens”*
929       | Same as above, except that two additional particles are
930         generated on the position of the oxygen, with names ‘LP1’ and
931         ‘LP2.’ This is for use with five-atom water models such as
932         TIP5P \ :ref:`130 <refMahoney2000a>`.
934 -  The name of the new H atom (or its prefix, *e.g.* ``HD2``
935    for the asparagine example given earlier).
937 -  Three or four control atoms (i,j,k,l), where the first always is the
938    atom to which the H atoms are connected. The other two or three
939    depend on the code selected. For water, there is only one control
940    atom.
942 Some more exotic cases can be approximately constructed from the above
943 tools, and with suitable use of energy minimization are good enough for
944 beginning MD simulations. For example secondary amine hydrogen, nitrenyl
945 hydrogen (:math:`\mathrm{C}=\mathrm{NH}`)
946 and even ethynyl hydrogen could be approximately constructed using
947 method 2 above for hydroxyl hydrogen.
949 Termini database
950 ~~~~~~~~~~~~~~~~
952 The termini
953 databases
954 are stored in ``aminoacids.n.tdb`` and
955 ``aminoacids.c.tdb`` for the N- and C-termini respectively.
956 They contain information for the :ref:`pdb2gmx <gmx pdb2gmx>` program on how
957 to connect new atoms to existing ones, which atoms should be removed or
958 changed, and which bonded interactions should be added. Their format is
959 as follows (from ``gromos43a1.ff/aminoacids.c.tdb``):
963     [ None ]
965     [ COO- ]
966     [ replace ]
967     C   C       C       12.011  0.27
968     O   O1      OM      15.9994 -0.635
969     OXT O2      OM      15.9994 -0.635
970     [ add ]
971     2   8       O       C       CA      N
972         OM      15.9994 -0.635
973     [ bonds ]
974     C   O1      gb_5
975     C   O2      gb_5
976     [ angles ]
977     O1  C       O2      ga_37
978     CA  C       O1      ga_21
979     CA  C       O2      ga_21
980     [ dihedrals ]
981     N   CA      C       O2      gd_20
982     [ impropers ]
983     C   CA      O2      O1      gi_1
985 The file is organized in blocks, each with a header specifying the name
986 of the block. These blocks correspond to different types of termini that
987 can be added to a molecule. In this example ``[ COO- ]`` is
988 the first block, corresponding to changing the terminal carbon atom into
989 a deprotonated carboxyl group. ``[ None ]`` is the second
990 terminus type, corresponding to a terminus that leaves the molecule as
991 it is. Block names cannot be any of the following:
992 ``replace``, ``add``, ``delete``,
993 ``bonds``, ``angles``,
994 ``dihedrals``, ``impropers``. Doing so would
995 interfere with the parameters of the block, and would probably also be
996 very confusing to human readers.
998 For each block the following options are present:
1000 -  | ``[ replace ]``
1001    | Replace an existing atom by one with a different atom type, atom
1002      name, charge, and/or mass. This entry can be used to replace an
1003      atom that is present both in the input coordinates and in the
1004      :ref:`rtp` database, but also to only rename an atom in
1005      the input coordinates such that it matches the name in the force
1006      field. In the latter case, there should also be a corresponding
1007      ``[ add ]`` section present that gives instructions to
1008      add the same atom, such that the position in the sequence and the
1009      bonding is known. Such an atom can be present in the input
1010      coordinates and kept, or not present and constructed by
1011      :ref:`pdb2gmx <gmx pdb2gmx>`. For each atom to be replaced on line
1012      should be entered with the following fields:
1014    -  name of the atom to be replaced
1016    -  new atom name (optional)
1018    -  new atom type
1020    -  new mass
1022    -  new charge
1024 -  | ``[ add ]``
1025    | Add new atoms. For each (group of) added atom(s), a two-line entry
1026      is necessary. The first line contains the same fields as an entry
1027      in the hydrogen database (name of the new atom, number of atoms,
1028      type of addition, control atoms, see :ref:`hdb`), but the
1029      possible types of addition are extended by two more, specifically
1030      for C-terminal additions:
1032    #. | *two carboxyl oxygens, -COO*:math:`^-`
1033       | Two oxygens (n1,n2) are generated according to rule 3, at a
1034         distance of 0.136 nm from atom i and an angle
1035         (n1-i-j)=(n2-i-j)=117 degrees
1037    #. | *carboxyl oxygens and hydrogen, -COOH*
1038       | Two oxygens (n1,n2) are generated according to rule 3, at
1039         distances of 0.123 nm and 0.125 nm from atom i for n1 and n2,
1040         respectively, and angles (n1-i-j)=121 and (n2-i-j)=115 degrees.
1041         One hydrogen (n:math:`^\prime`) is generated around n2 according
1042         to rule 2, where n-i-j and n-i-j-k should be read as
1043         n\ :math:`^\prime`-n2-i and n\ :math:`^\prime`-n2-i-j,
1044         respectively.
1046    After this line, another line follows that specifies the details of
1047    the added atom(s), in the same way as for replacing atoms, *i.e.*:
1049    -  atom type
1051    -  mass
1053    -  charge
1055    -  charge group (optional)
1057    Like in the hydrogen database (see :ref:`rtp`), when more than one
1058    atom is connected to an existing one, a number will be appended to
1059    the end of the atom name. **Note** that, like in the hydrogen
1060    database, the atom name is now on the same line as the control atoms,
1061    whereas it was at the beginning of the second line prior to |Gromacs|
1062    version 3.3. When the charge group field is left out, the added atom
1063    will have the same charge group number as the atom that it is bonded
1064    to.
1066 -  | ``[ delete ]``
1067    | Delete existing atoms. One atom name per line.
1069 -  | ``[ bonds ]``, ``[ angles ]``,
1070      ``[ dihedrals ]`` and ``[ impropers ]``
1071    | Add additional bonded parameters. The format is identical to that
1072      used in the :ref:`rtp` file, see :ref:`rtp`.
1074 Virtual site database
1075 ~~~~~~~~~~~~~~~~~~~~~
1077 Since we cannot rely on the positions of hydrogens in input files, we
1078 need a special input file to decide the geometries and parameters with
1079 which to add virtual site hydrogens. For more complex virtual site
1080 constructs (*e.g.* when entire aromatic side chains are made rigid) we
1081 also need information about the equilibrium bond lengths and angles for
1082 all atoms in the side chain. This information is specified in the
1083 :ref:`vsd` file for each force field. Just as for the termini,
1084 there is one such file for each class of residues in the
1085 :ref:`rtp` file.
1087 The virtual site database is not really a very simple list of
1088 information. The first couple of sections specify which mass centers
1089 (typically called MCH\ :math:`_3`/MNH:math:`_3`) to use for
1090 CH\ :math:`_3`, NH\ :math:`_3`, and NH\ :math:`_2` groups. Depending on
1091 the equilibrium bond lengths and angles between the hydrogens and heavy
1092 atoms we need to apply slightly different constraint distances between
1093 these mass centers. **Note** that we do *not* have to specify the actual
1094 parameters (that is automatic), just the type of mass center to use. To
1095 accomplish this, there are three sections names ``[ CH3 ]``,
1096 ``[ NH3 ]``, and ``[ NH2 ]``. For each of these we expect three columns.
1097 The first column is the atom type bound to the 2/3 hydrogens, the second
1098 column is the next heavy atom type which this is bound, and the third
1099 column the type of mass center to use. As a special case, in the
1100 ``[ NH2 ]`` section it is also possible to specify ``planar`` in the
1101 second column, which will use a different construction without mass
1102 center. There are currently different opinions in some force fields
1103 whether an NH\ :math:`_2` group should be planar or not, but we try hard
1104 to stick to the default equilibrium parameters of the force field.
1106 The second part of the virtual site database contains explicit
1107 equilibrium bond lengths and angles for pairs/triplets of atoms in
1108 aromatic side chains. These entries are currently read by specific
1109 routines in the virtual site generation code, so if you would like to
1110 extend it *e.g.* to nucleic acids you would also need to write new code
1111 there. These sections are named after the short amino acid names
1112 (``[ PHE ]``, ``[ TYR ]``, ``[ TRP ]``, ``[ HID ]``, ``[ HIE ]``,
1113 ``[ HIP ]``), and simply contain 2 or 3 columns with atom names,
1114 followed by a number specifying the bond length (in nm) or angle (in
1115 degrees). **Note** that these are approximations of the equilibrated
1116 geometry for the entire molecule, which might not be identical to the
1117 equilibrium value for a single bond/angle if the molecule is strained.
1119 .. _specbond:
1121 Special bonds
1122 ~~~~~~~~~~~~~
1124 The primary mechanism used by
1125 :ref:`pdb2gmx <gmx pdb2gmx>` to generate
1126 inter-residue bonds relies on head-to-tail linking of backbone atoms in
1127 different residues to build a macromolecule. In some cases (*e.g.*
1128 disulfide bonds, a heme
1129 group, branched
1130 polymers), it is necessary to
1131 create inter-residue bonds that do not lie on the backbone. The file
1132 ``specbond.dat`` takes
1133 care of this function. It is necessary that the residues belong to the
1134 same ``[ moleculetype ]``. The ``-merge`` and
1135 ``-chainsep`` functions of :ref:`pdb2gmx <gmx pdb2gmx>` can be
1136 useful when managing special inter-residue bonds between different
1137 chains.
1139 The first line of ``specbond.dat`` indicates the number of
1140 entries that are in the file. If you add a new entry, be sure to
1141 increment this number. The remaining lines in the file provide the
1142 specifications for creating bonds. The format of the lines is as
1143 follows:
1145 ``resA atomA nbondsA resB atomB nbondsB length newresA
1146 newresB``
1148 The columns indicate:
1150 #. ``resA`` The name of residue A that participates in the
1151    bond.
1153 #. ``atomA`` The name of the atom in residue A that forms
1154    the bond.
1156 #. ``nbondsA`` The total number of bonds
1157    ``atomA`` can form.
1159 #. ``resB`` The name of residue B that participates in the
1160    bond.
1162 #. ``atomB`` The name of the atom in residue B that forms
1163    the bond.
1165 #. ``nbondsB`` The total number of bonds
1166    ``atomB`` can form.
1168 #. ``length`` The reference length for the bond. If
1169    ``atomA`` and ``atomB`` are not within
1170    ``length`` :math:`\pm` 10% in the coordinate file
1171    supplied to :ref:`pdb2gmx <gmx pdb2gmx>`, no bond will be formed.
1173 #. ``newresA`` The new name of residue A, if necessary. Some
1174    force fields use *e.g.* CYS2 for a cysteine in a disulfide or heme
1175    linkage.
1177 #. ``newresB`` The new name of residue B, likewise.
1179 File formats
1180 ------------
1182 .. _topfile:
1184 Topology file
1185 ~~~~~~~~~~~~~
1187 The topology file is built following the |Gromacs| specification for a
1188 molecular topology. A :ref:`top` file can be generated by
1189 :ref:`pdb2gmx <gmx pdb2gmx>`. All possible entries in the topology file are
1190 listed in :numref:`Tables %s <tab-topfile1>` and
1191 :numref:`%s <tab-topfile2>`. Also tabulated are: all the units of
1192 the parameters, which interactions can be perturbed for free energy
1193 calculations, which bonded interactions are used by
1194 :ref:`grompp <gmx grompp>` for generating exclusions, and which bonded
1195 interactions can be converted to constraints by :ref:`grompp <gmx grompp>`.
1197 .. |VCR| replace:: V\ :math:`^{(cr)}`
1198 .. |WCR| replace:: W\ :math:`^{(cr)}`
1199 .. |CRO| replace:: :math:`^{(cr)}`
1200 .. |TREF| replace:: :numref:`Table %s <tab-topfile2>`
1201 .. |AKJM| replace:: :math:`a~\mathrm{kJ~mol}^{-1}`
1202 .. |KJN6| replace:: :math:`\mathrm{kJ~mol}^{-1}~\mathrm{nm}^{-6}`
1203 .. |BNM| replace:: :math:`b~\mathrm{nm}^{-1}`
1204 .. |C6LJ| replace:: :math:`c_6`
1205 .. |STAR| replace:: :math:`^{(*)}`
1206 .. |NREX| replace:: :math:`n_{ex}^{(nrexcl)}`
1207 .. |QEMU| replace:: :math:`q` (e); :math:`m` (u) 
1208 .. |MQM| replace:: :math:`q,m`
1210 .. _tab-topfile1:
1212 .. table:: The :ref:`topology <top>` file.
1214         +------------------------------------------------------------------------------------------------------------+
1215         | Parameters                                                                                                 |
1216         +===================+===========================+=====+====+=========================================+=======+
1217         | interaction type  | directive                 | #   | f. | parameters                              | F. E. |
1218         |                   |                           | at. | tp |                                         |       |
1219         +-------------------+---------------------------+-----+----+-----------------------------------------+-------+
1220         | *mandatory*       | ``defaults``              |            non-bonded function type;                       |
1221         |                   |                           |            combination rule\ |CRO|;                        |
1222         |                   |                           |            generate pairs (no/yes);                        |
1223         |                   |                           |            fudge LJ (); fudge QQ ()                        |
1224         +-------------------+---------------------------+------------------------------------------------------------+
1225         | *mandatory*       | ``atomtypes``             |            atom type; m (u); q (e); particle type;         | 
1226         |                   |                           |            |VCR| ; |WCR|                                   |
1227         +-------------------+---------------------------+------------------------------------------------------------+
1228         |                   | ``bondtypes``             |  (see |TREF|, directive ``bonds``)                         |
1229         +                   +                           +                                                            +
1230         |                   | ``pairtypes``             |  (see |TREF|, directive ``pairs``)                         |
1231         +                   +                           +                                                            +
1232         |                   | ``angletypes``            |  (see |TREF|, directive ``angles``)                        |
1233         +                   +                           +                                                            +
1234         |                   | ``dihedraltypes``\ |STAR| |  (see |TREF|, directive ``dihedrals``)                     |
1235         +                   +                           +                                                            +
1236         |                   | ``constrainttypes``       |  (see |TREF|, directive ``constraints``)                   |
1237         +-------------------+---------------------------+-----+----+-------------------------------------------------+
1238         | LJ                | ``nonbond_params``        |  2  | 1  |  |VCR|  ; |WCR|                                 |
1239         +                   +                           +     +    +                                                 +
1240         | Buckingham        | ``nonbond_params``        |  2  | 2  |  |AKJM| ; |BNM|;                                |
1241         |                   |                           |     |    |  |C6LJ| (|KJN6|)                                |
1242         +-------------------+---------------------------+-----+----+-------------------------------------------------+
1244 .. table:: 
1246         +------------------------------------------------------------------------------------------------------------+
1247         | Molecule definition(s)                                                                                     |
1248         +===================+===========================+============================================================+
1249         | *mandatory*       | ``moleculetype``          | molecule name; |NREX|                                      |
1250         +-------------------+---------------------------+-----+----------------------------------------------+-------+
1251         | *mandatory*       | ``atoms``                 | 1   | atom type; residue number;                   | type  |
1252         |                   |                           |     | residue name; atom name;                     |       |
1253         |                   |                           |     | charge group number; |QEMU|                  | |MQM| |
1254         +-------------------+---------------------------+-----+----------------------------------------------+-------+
1255         | intra-molecular interaction and geometry definitions as described in |TREF|                                |
1256         +------------------------------------------------------------------------------------------------------------+
1258 .. table::
1260         +-------------+---------------+------------------------------------+
1261         | System      |               |                                    |
1262         +=============+===============+====================================+
1263         | *mandatory* | ``system``    | system name                        |
1264         +-------------+---------------+------------------------------------+
1265         | *mandatory* | ``molecules`` | molecule name; number of molecules |
1266         +-------------+---------------+------------------------------------+
1268 .. table::
1270         +------------------------------+----------------------------------------------------+
1271         | Inter-molecular interactions |                                                    |
1272         +==============================+====================================================+
1273         | *optional*                   | ``intermolecular_interactions``                    |
1274         +------------------------------+----------------------------------------------------+
1275         | one or more bonded interactions as described in |TREF|, with two or more atoms,   |
1276         | no interactions that generate exclusions, no constraints, use global atom numbers |
1277         +-----------------------------------------------------------------------------------+
1279 .. parsed-literal::
1281     '\# at' is the required number of atom type indices for this directive
1282     'f. tp' is the value used to select this function type
1283     'F. E.' indicates which of the parameters can be interpolated in free energy calculations
1284     |CRO| the combination rule determines the type of LJ parameters, see 
1285     |STAR| for ``dihedraltypes`` one can specify 4 atoms or the inner (outer for improper) 2 atoms
1286     |NREX| exclude neighbors :math:`n_{ex}` bonds away for non-bonded interactions
1287     For free energy calculations, type, :math:`q` and :math:`m`  or no parameters should be added
1288     for topology 'B' (:math:`\lambda = 1`) on the same line, after the normal parameters.
1290 .. |BZERO| replace:: :math:`b_0`
1291 .. |KB| replace:: :math:`k_b`
1292 .. |KDR| replace:: :math:`k_{dr}`
1293 .. |NM2| replace:: (kJ mol\ :math:`^{-1}`\ nm\ :math:`^{-2}`
1294 .. |NM4| replace:: (kJ mol\ :math:`^{-1}`\ nm\ :math:`^{-4}`
1295 .. |DKJ| replace:: :math:`D` (kJ mol\ :math:`^{-1}`
1296 .. |BETA| replace:: :math:`\beta` (nm\ :math:`^{-1}`
1297 .. |C23| replace:: :math:`C_{i=2,3}` (kJ mol\ :math:`^{-1}\ nm\ :math:`^{-i}`
1298 .. |BMM| replace:: :math:`b_m`
1299 .. |GE0| replace:: :math:`\geq 0`
1300 .. |KO| replace:: :math:`k` 
1301 .. |KJM| replace:: kJ mol\ :math:`^{-1}`
1302 .. |LUU| replace:: low, up\ :math:`_1`,\ :math:`_2`
1303 .. |MV| replace:: :math:`V`
1304 .. |MW| replace:: :math:`W`
1305 .. |QIJ| replace:: :math:`q_i`; :math:`q_j`
1306 .. |THE0| replace:: :math:`\theta_0`
1307 .. |KTHE| replace:: :math:`k_\theta`
1308 .. |KJR2| replace:: kJ mol\ :math:`^{-1}`\ rad\ :math:`^{-2}`
1309 .. |RN13| replace:: :math:`r_{13}`
1310 .. |KUB| replace:: :math:`k_{UB}`
1311 .. |C024| replace:: :math:`C_{i=0,1,2,3,4}`
1312 .. |KJRI| replace:: kJ mol\ :math:`^{-1}`\ rad\ :math:`^{-i}`
1313 .. |PHIS| replace:: :math:`\phi_s`
1314 .. |PHI0| replace:: :math:`\phi_0`
1315 .. |KPHI| replace:: :math:`k_\phi`
1316 .. |PHIK| replace:: :math:`\phi,k`
1317 .. |XI0| replace:: :math:`\xi_0`
1318 .. |KXI| replace:: :math:`k_\xi`
1319 .. |C0| replace:: :math:`C_0`
1320 .. |C1| replace:: :math:`C_1`
1321 .. |C2| replace:: :math:`C_2`
1322 .. |C3| replace:: :math:`C_3`
1323 .. |C4| replace:: :math:`C_4`
1324 .. |C5| replace:: :math:`C_5`
1325 .. |A0| replace:: :math:`a_0`
1326 .. |A1| replace:: :math:`a_1`
1327 .. |A2| replace:: :math:`a_2`
1328 .. |A3| replace:: :math:`a_3`
1329 .. |A4| replace:: :math:`a_4`
1330 .. |DOH| replace:: :math:d_{\mbox{\sc oh}}`
1331 .. |DHH| replace:: :math:d_{\mbox{\sc hh}}`
1332 .. |AO| replace:: :math:`a`
1333 .. |BO| replace:: :math:`b`
1334 .. |CO| replace:: :math:`c`
1335 .. |DO| replace:: :math:`d`
1336 .. |KX| replace:: :math:`k_{x}`
1337 .. |KY| replace:: :math:`k_{y}`
1338 .. |KZ| replace:: :math:`k_{z}`
1339 .. |GO| replace:: :math:`g`
1340 .. |RO| replace:: :math:`r`
1341 .. |DPHI| replace:: :math:`\Delta\phi`
1342 .. |DIHR| replace:: :math:`k_{\mathrm{dihr}}`
1343 .. |THET| replace:: :math:`\theta`
1344 .. |NM| replace:: nm\ :math:`^{-1}`
1345 .. |KC| replace:: :math:`k_c`
1346 .. |THEK| replace:: :math:`\theta,k`
1347 .. |R1E| replace:: :math:`r_{1e}`
1348 .. |R2E| replace:: :math:`r_{2e}`
1349 .. |R3E| replace:: :math:`r_{3e}`
1350 .. |KRR| replace:: :math:`k_{rr'}`
1351 .. |KRTH| replace:: :math:`k_{r\theta}`
1352 .. |ALPH| replace:: :math:`\alpha`; |CO| (U nm\ :math:`^{\alpha}`
1353 .. |UM1| replace:: U\ :math:`^{-1}`
1355 .. _tab-topfile2:
1357 .. table:: Details of ``[ moleculetype ]`` directives
1359             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1360             | Name of interaction                | Topology file directive    | num.       | func.     | Order of parameters and their units                                     | use in     | 
1361             |                                    |                            | atoms [1]_ | type [2]_ |                                                                         | F.E.? [3]_ |
1362             +====================================+============================+============+===========+=========================================================================+============+
1363             | bond                               | ``bonds`` [4]_, [5]_       | 2          | 1         | |BZERO| (nm); |KB| |NM2|                                                | all        | 
1364             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1365             | G96 bond                           | ``bonds`` [4]_, [5]_       | 2          | 2         | |BZERO| (nm); |KB| |NM4|                                                | all        |
1366             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1367             | Morse                              | ``bonds`` [4]_, [5]_       | 2          | 3         | |BZERO| (nm); |DKJ|; |BETA|                                             | all        | 
1368             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1369             | cubic bond                         | ``bonds`` [4]_, [5]_       | 2          | 4         | |BZERO| (nm); |C23|                                                     |            | 
1370             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1371             | connection                         | ``bonds`` [4]_             | 2          | 5         |                                                                         |            | 
1372             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1373             | harmonic potential                 | ``bonds``                  | 2          | 6         | |BZERO| (nm); |KB| |NM2|                                                | all        | 
1374             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1375             | FENE bond                          | ``bonds`` [4]_             | 2          | 7         | |BMM|   (nm); |KB| |NM2|                                                |            | 
1376             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1377             | tabulated bond                     | ``bonds`` [4]_             | 2          | 8         | table number (|GE0|); |KO| |KJM|                                        | |KO|       |
1378             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1379             | tabulated bond [6]_                | ``bonds``                  | 2          | 9         | table number (|GE0|); |KO| |KJM|                                        | |KO|       |
1380             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1381             | restraint potential                | ``bonds``                  | 2          | 10        | |LUU| (nm); |KDR| (|NM2|)                                               | all        | 
1382             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1383             | extra LJ or Coulomb                | ``pairs``                  | 2          | 1         | |MV| [7]_; |MW| [7]_                                                    | all        | 
1384             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1385             | extra LJ or Coulomb                | ``pairs``                  | 2          | 2         | fudge QQ (); |QIJ| (e), |MV| [7]_; |MW| [7]_                            |            | 
1386             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1387             | extra LJ or Coulomb                | ``pairs_nb``               | 2          | 1         | |QIJ| (e); |MV| [7]_; |MW| [7]_                                         |            | 
1388             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1389             | angle                              | ``angles`` [5]_            | 3          | 1         | |THE0| (deg); |KTHE| (|KJR2|)                                           | all        | 
1390             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1391             | G96 angle                          | ``angles`` [5]_            | 3          | 2         | |THE0| (deg); |KTHE| (|KJM|)                                            | all        | 
1392             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1393             | cross bond-bond                    | ``angles``                 | 3          | 3         | |R1E|, |R2E| (nm); |KRR| (|NM2|)                                        |            | 
1394             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1395             | cross bond-angle                   | ``angles``                 | 3          | 4         | |R1E|, |R2E|, |R3E| (nm); |KRTH| (|NM2|)                                |            | 
1396             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1397             | Urey-Bradley                       | ``angles`` [5]_            | 3          | 5         | |THE0| (deg); |KTHE| (|KJR2|); |RN13| (nm); |KUB| (|NM2|)               | all        |
1398             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1399             | quartic angle                      | ``angles`` [5]_            | 3          | 6         | |THE0| (deg); |C024| (|KJRI|)                                           |            | 
1400             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1401             | tabulated angle                    | ``angles``                 | 3          | 8         | table number (|GE0|); |KO| (|KJM|)                                      | |KO|       | 
1402             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1403             |  |  restricted                     |                            |            |           |                                                                         |            |
1404             |  |  bending potential              | ``angles``                 | 3          | 10        | |THE0| (deg); |KTHE| (|KJM|)                                            |            | 
1405             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1406             | proper dihedral                    | ``dihedrals``              | 4          | 1         | |PHIS| (deg); |KPHI| (|KJM|); multiplicity                              | |PHIK|     | 
1407             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1408             | improper dihedral                  | ``dihedrals``              | 4          | 2         | |XI0| (deg); |KXI| (|KJR2|)                                             | all        | 
1409             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1410             | Ryckaert-Bellemans dihedral        | ``dihedrals``              | 4          | 3         | |C0|, |C1|, |C2|, |C3|, |C4|, |C5| (|KJM|)                              | all        | 
1411             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1412             | periodic improper dihedral         | ``dihedrals``              | 4          | 4         | |PHIS| (deg); |KPHI| (|KJM|); multiplicity                              | |PHIK|     | 
1413             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1414             | Fourier dihedral                   | ``dihedrals``              | 4          | 5         | |C1|, |C2|, |C3|, |C4|, |C5| (|KJM|)                                    | all        | 
1415             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1416             | tabulated dihedral                 | ``dihedrals``              | 4          | 8         | table number (|GE0|); |KO| (|KJM|)                                      | |KO|       |
1417             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1418             | proper dihedral (multiple)         | ``dihedrals``              | 4          | 9         | |PHIS| (deg); |KPHI| (|KJM|); multiplicity                              | |PHIK|     | 
1419             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1420             | restricted dihedral                | ``dihedrals``              | 4          | 10        | |PHI0| (deg); |KPHI| (|KJM|)                                            |            | 
1421             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1422             | combined bending-torsion potential | ``dihedrals``              | 4          | 11        | |A0|, |A1|, |A2|, |A3|, |A4| (|KJM|)                                    |            | 
1423             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1424             | exclusions                         | ``exclusions``             | 1          |           | one or more atom indices                                                |            | 
1425             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1426             | constraint                         | ``constraints`` [4]_       | 2          | 1         | |BZERO| (nm)                                                            | all        | 
1427             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1428             | constraint [6]_                    | ``constraints``            | 2          | 2         | |BZERO| (nm)                                                            | all        | 
1429             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1430             | SETTLE                             | ``settles``                | 1          | 1         | |DOH|, |DHH| (nm)                                                       |            | 
1431             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1432             | 2-body virtual site                | ``virtual_sites2``         | 3          | 1         | |AO| ()                                                                 |            | 
1433             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1434             | 3-body virtual site                | ``virtual_sites3``         | 4          | 1         | |AO|, |BO| ()                                                           |            | 
1435             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1436             | 3-body virtual site (fd)           | ``virtual_sites3``         | 4          | 2         | |AO| (); |DO| (nm)                                                      |            | 
1437             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1438             | 3-body virtual site (fad)          | ``virtual_sites3``         | 4          | 3         | |THET| (deg); |DO| (nm)                                                 |            | 
1439             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1440             | 3-body virtual site (out)          | ``virtual_sites3``         | 4          | 4         | |AO|, |BO| (); |CO| (|NM|)                                              |            | 
1441             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1442             | 4-body virtual site (fdn)          | ``virtual_sites4``         | 5          | 2         | |AO|, |BO| (); |CO| (nm)                                                |            | 
1443             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1444             | N-body virtual site (COG)          | ``virtual_sitesn``         | 1          | 1         | one or more constructing atom indices                                   |            | 
1445             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1446             | N-body virtual site (COM)          | ``virtual_sitesn``         | 1          | 2         | one or more constructing atom indices                                   |            | 
1447             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1448             | N-body virtual site (COW)          | ``virtual_sitesn``         | 1          | 3         |  |  one or more pairs consisting of                                     |            |
1449             |                                    |                            |            |           |  |  constructing atom index and weight                                  |            | 
1450             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1451             | position restraint                 | ``position_restraints``    | 1          | 1         | |KX|, |KY|, |KZ| (|NM2|)                                                | all        |
1452             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1453             | flat-bottomed position restraint   | ``position_restraints``    | 1          | 2         | |GO|, |RO| (nm), |KO| (|NM2|)                                           |            | 
1454             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1455             | distance restraint                 | ``distance_restraints``    | 2          | 1         | type; label; |LUU| (nm); weight ()                                      |            | 
1456             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1457             | dihedral restraint                 | ``dihedral_restraints``    | 4          | 1         | |PHI0| (deg); |DPHI| (deg); |DIHR| (|KJR2|)                             | all        | 
1458             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1459             | orientation restraint              | ``orientation_restraints`` | 2          | 1         | exp.; label; |ALPH|; obs. (U); weight (|UM1|)                           |            |
1460             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1461             | angle restraint                    | ``angle_restraints``       | 4          | 1         | |THE0| (deg); |KC| (|KJM|); multiplicity                                | |THEK|     | 
1462             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1463             | angle restraint (z)                | ``angle_restraints_z``     | 2          | 1         | |THE0| (deg); |KC| (|KJM|); multiplicity                                | |THEK|     | 
1464             +------------------------------------+----------------------------+------------+-----------+-------------------------------------------------------------------------+------------+
1466 .. [1]
1467    The required number of atom indices for this directive
1468    
1469 .. [2]
1470    The index to use to select this function type
1471    
1472 .. [3]
1473    Indicates which of the parameters can be interpolated in free energy calculations
1474    
1475 .. [4]
1476    This interaction type will be used by :ref:`grompp <gmx grompp>` for generating exclusions
1477    
1478 .. [5]
1479    This interaction type can be converted to constraints by :ref:`grompp <gmx grompp>`
1480    
1481 .. [7]
1482    The combination rule determines the type of LJ parameters, see
1483    
1484 .. [6]
1485    No connection, and so no exclusions, are generated for this interaction
1487 Description of the file layout:
1489 -  Semicolon (;) and newline characters surround comments
1491 -  On a line ending with :math:`\backslash` the newline character is
1492    ignored.
1494 -  Directives are surrounded by ``[`` and ``]``
1496 -  The topology hierarchy (which must be followed) consists of three
1497    levels:
1499    -  the parameter level, which defines certain force-field
1500       specifications (see :numref:`Table %s <tab-topfile1>`)
1502    -  the molecule level, which should contain one or more molecule
1503       definitions (see :numref:`Table %s <tab-topfile2>`)
1505    -  the system level, containing only system-specific information
1506       (``[ system ]`` and ``[ molecules ]``)
1508 -  Items should be separated by spaces or tabs, not commas
1510 -  Atoms in molecules should be numbered consecutively starting at 1
1512 -  Atoms in the same charge group must be listed consecutively
1514 -  The file is parsed only once, which implies that no forward
1515    references can be treated: items must be defined before they can be
1516    used
1518 -  Exclusions can be generated from the bonds or overridden manually
1520 -  The bonded force types can be generated from the atom types or
1521    overridden per bond
1523 -  It is possible to apply multiple bonded interactions of the same type
1524    on the same atoms
1526 -  Descriptive comment lines and empty lines are highly recommended
1528 -  Starting with |Gromacs| version 3.1.3, all directives at the parameter
1529    level can be used multiple times and there are no restrictions on the
1530    order, except that an atom type needs to be defined before it can be
1531    used in other parameter definitions
1533 -  If parameters for a certain interaction are defined multiple times
1534    for the same combination of atom types the last definition is used;
1535    starting with |Gromacs| version 3.1.3 :ref:`grompp <gmx grompp>` generates
1536    a warning for parameter redefinitions with different values
1538 -  Using one of the ``[ atoms ]``,
1539    ``[ bonds ]``, ``[ pairs ]``,
1540    ``[ angles ]``, etc. without having used
1541    ``[ moleculetype ]`` before is meaningless and generates
1542    a warning
1544 -  Using ``[ molecules ]`` without having used
1545    ``[ system ]`` before is meaningless and generates a
1546    warning.
1548 -  After ``[ system ]`` the only allowed directive is
1549    ``[ molecules ]``
1551 -  Using an unknown string in ``[ ]`` causes all the data
1552    until the next directive to be ignored and generates a warning
1554 Here is an example of a topology file, ``urea.top``:
1558     ;
1559     ;       Example topology file
1560     ;
1561     ; The force-field files to be included
1562     #include "amber99.ff/forcefield.itp"
1564     [ moleculetype ]
1565     ; name  nrexcl
1566     Urea         3
1568     [ atoms ]
1569        1  C  1  URE      C      1     0.880229  12.01000   ; amber C  type
1570        2  O  1  URE      O      2    -0.613359  16.00000   ; amber O  type
1571        3  N  1  URE     N1      3    -0.923545  14.01000   ; amber N  type
1572        4  H  1  URE    H11      4     0.395055   1.00800   ; amber H  type
1573        5  H  1  URE    H12      5     0.395055   1.00800   ; amber H  type
1574        6  N  1  URE     N2      6    -0.923545  14.01000   ; amber N  type
1575        7  H  1  URE    H21      7     0.395055   1.00800   ; amber H  type
1576        8  H  1  URE    H22      8     0.395055   1.00800   ; amber H  type
1578     [ bonds ]
1579         1       2
1580         1       3       
1581         1   6
1582         3       4
1583         3       5
1584         6       7
1585         6       8
1587     [ dihedrals ] 
1588     ;   ai    aj    ak    al funct  definition
1589          2     1     3     4   9     
1590          2     1     3     5   9     
1591          2     1     6     7   9     
1592          2     1     6     8   9     
1593          3     1     6     7   9     
1594          3     1     6     8   9     
1595          6     1     3     4   9     
1596          6     1     3     5   9     
1598     [ dihedrals ] 
1599          3     6     1     2   4     
1600          1     4     3     5   4         
1601          1     7     6     8   4
1603     [ position_restraints ]
1604     ; you wouldn't normally use this for a molecule like Urea,
1605     ; but we include it here for didactic purposes
1606     ; ai   funct    fc
1607        1     1     1000    1000    1000 ; Restrain to a point
1608        2     1     1000       0    1000 ; Restrain to a line (Y-axis)
1609        3     1     1000       0       0 ; Restrain to a plane (Y-Z-plane)
1611     [ dihedral_restraints ]
1612     ; ai   aj    ak    al  type  phi  dphi  fc
1613         3    6     1    2     1  180     0  10
1614         1    4     3    5     1  180     0  10
1616     ; Include TIP3P water topology
1617     #include "amber99/tip3p.itp"
1619     [ system ]
1620     Urea in Water
1622     [ molecules ]
1623     ;molecule name   nr.
1624     Urea             1
1625     SOL              1000
1627 Here follows the explanatory text.
1629 **#include “amber99.ff/forcefield.itp” :** this includes
1630 the information for the force field you are using, including bonded and
1631 non-bonded parameters. This example uses the AMBER99 force field, but
1632 your simulation may use a different force field. :ref:`grompp <gmx grompp>`
1633 will automatically go and find this file and copy-and-paste its content.
1634 That content can be seen in
1635 ``share/top/amber99.ff/forcefield.itp}``, and it
1640     #define _FF_AMBER
1641     #define _FF_AMBER99
1643     [ defaults ]
1644     ; nbfunc        comb-rule       gen-pairs       fudgeLJ fudgeQQ
1645     1               2               yes             0.5     0.8333
1647     #include "ffnonbonded.itp"
1648     #include "ffbonded.itp"
1650 The two ``#define`` statements set up the conditions so that
1651 future parts of the topology can know that the AMBER 99 force field is
1652 in use.
1654 **[ defaults ] :**
1656 -  ``nbfunc`` is the non-bonded function type. Use 1 (Lennard-Jones) or 2
1657    (Buckingham)
1659 -  ``comb-rule`` is the number of the combination rule (see :ref:`nbpar`).
1661 -  ``gen-pairs`` is for pair generation. The default is
1662    ‘no’, *i.e.* get 1-4 parameters from the pairtypes list. When
1663    parameters are not present in the list, stop with a fatal error.
1664    Setting ‘yes’ generates 1-4 parameters that are not present in the
1665    pair list from normal Lennard-Jones parameters using
1666    ``fudgeLJ``
1668 -  ``fudgeLJ`` is the factor by which to multiply
1669    Lennard-Jones 1-4 interactions, default 1
1671 -  ``fudgeQQ`` is the factor by which to multiply
1672    electrostatic 1-4 interactions, default 1
1674 -  :math:`N` is the power for the repulsion term in a 6-\ :math:`N`
1675    potential (with nonbonded-type Lennard-Jones only), starting with
1676    |Gromacs| version 4.5, :ref:`grompp <gmx mdrun>` also reads and applies
1677    :math:`N`, for values not equal to 12 tabulated interaction functions
1678    are used (in older version you would have to use user tabulated
1679    interactions).
1681 **Note** that ``gen-pairs``, ``fudgeLJ``,
1682 ``fudgeQQ``, and :math:`N` are optional.
1683 ``fudgeLJ`` is only used when generate pairs is set to
1684 ‘yes’, and ``fudgeQQ`` is always used. However, if you want
1685 to specify :math:`N` you need to give a value for the other parameters
1686 as well.
1688 Then some other ``#include`` statements add in the large
1689 amount of data needed to describe the rest of the force field. We will
1690 skip these and return to ``urea.top``. There we will see
1692 **[ moleculetype ] :** defines the name of your molecule
1693 in this :ref:`top` and nrexcl = 3 stands for excluding
1694 non-bonded interactions between atoms that are no further than 3 bonds
1695 away.
1697 **[ atoms ] :** defines the molecule, where
1698 ``nr`` and ``type`` are fixed, the rest is user
1699 defined. So ``atom`` can be named as you like,
1700 ``cgnr`` made larger or smaller (if possible, the total
1701 charge of a charge group should be zero), and charges can be changed
1702 here too.
1704 **[ bonds ] :** no comment.
1706 **[ pairs ] :** LJ and Coulomb 1-4 interactions
1708 **[ angles ] :** no comment
1710 **[ dihedrals ] :** in this case there are 9 proper
1711 dihedrals (funct = 1), 3 improper (funct = 4) and no Ryckaert-Bellemans
1712 type dihedrals. If you want to include Ryckaert-Bellemans type dihedrals
1713 in a topology, do the following (in case of *e.g.* decane):
1717     [ dihedrals ]
1718     ;  ai    aj    ak    al funct       c0       c1       c2
1719         1    2     3     4     3 
1720         2    3     4     5     3
1722 In the original implementation of the potential for
1723 alkanes \ :ref:`131 <refRyckaert78>` no 1-4 interactions were used, which means that in
1724 order to implement that particular force field you need to remove the
1725 1-4 interactions from the ``[ pairs ]`` section of your
1726 topology. In most modern force fields, like OPLS/AA or Amber the rules
1727 are different, and the Ryckaert-Bellemans potential is used as a cosine
1728 series in combination with 1-4 interactions.
1730 **[ position_restraints ] :** harmonically restrain the selected particles to reference
1731 positions (:ref:`positionrestraint`). The reference positions are read
1732 from a separate coordinate file by :ref:`grompp <gmx grompp>`.
1734 **[ dihedral_restraints ] :** restrain selected dihedrals to a reference value. The
1735 implementation of dihedral restraints is described in section
1736 :ref:`dihedralrestraint` of the manual. The parameters specified in
1737 the ``[dihedral_restraints]`` directive are as follows:
1739 -  ``type`` has only one possible value which is 1
1741 -  ``phi`` is the value of :math:`\phi_0` in :eq:`eqn. %s <eqndphi>` and
1742    :eq:`eqn. %s <eqndihre>` of the manual.
1744 -  ``dphi`` is the value of :math:`\Delta\phi` in :eq:`eqn. %s <eqndihre>` of the
1745    manual.
1747 -  ``fc`` is the force constant :math:`k_{dihr}` in :eq:`eqn. %s <eqndihre>` of the
1748    manual.
1750 **#include “tip3p.itp” :** includes a topology file that was already
1751 constructed (see section :ref:`molitp`).
1753 **[ system ] :** title of your system, user-defined
1755 **[ molecules ] :** this defines the total number of (sub)molecules in your system
1756 that are defined in this :ref:`top`. In this example file, it stands for 1
1757 urea molecule dissolved in 1000 water molecules. The molecule type ``SOL``
1758 is defined in the ``tip3p.itp`` file. Each name here must correspond to a
1759 name given with ``[ moleculetype ]`` earlier in the topology. The order of the blocks of
1760 molecule types and the numbers of such molecules must match the
1761 coordinate file that accompanies the topology when supplied to :ref:`grompp <gmx grompp>`.
1762 The blocks of molecules do not need to be contiguous, but some tools
1763 (e.g. :ref:`genion <gmx genion>`) may act only on the first or last such block of a
1764 particular molecule type. Also, these blocks have nothing to do with the
1765 definition of groups (see sec. :ref:`groupconcept` and
1766 sec. :ref:`usinggroups`).
1768 .. _molitp:
1770 Molecule.itp file
1771 ~~~~~~~~~~~~~~~~~
1773 If you construct a topology file you will use frequently (like the water
1774 molecule, ``tip3p.itp``, which is already constructed for
1775 you) it is good to make a ``molecule.itp`` file. This only
1776 lists the information of one particular molecule and allows you to
1777 re-use the ``[ moleculetype ]`` in multiple systems without
1778 re-invoking :ref:`pdb2gmx <gmx pdb2gmx>` or manually copying and pasting. An
1779 example ``urea.itp`` follows:
1783     [ moleculetype ]
1784     ; molname   nrexcl
1785     URE         3
1787     [ atoms ]
1788        1  C  1  URE      C      1     0.880229  12.01000   ; amber C  type
1789     ...
1790        8  H  1  URE    H22      8     0.395055   1.00800   ; amber H  type
1792     [ bonds ]
1793         1       2
1794     ...
1795         6       8
1796     [ dihedrals ] 
1797     ;   ai    aj    ak    al funct  definition
1798          2     1     3     4   9     
1799     ...
1800          6     1     3     5   9     
1801     [ dihedrals ] 
1802          3     6     1     2   4     
1803          1     4     3     5   4         
1804          1     7     6     8   4
1806 Using :ref:`itp` files results in a very short
1807 :ref:`top` file:
1811     ;
1812     ;       Example topology file
1813     ;
1814     ; The force field files to be included
1815     #include "amber99.ff/forcefield.itp"
1817     #include "urea.itp"
1819     ; Include TIP3P water topology
1820     #include "amber99/tip3p.itp"
1822     [ system ]
1823     Urea in Water
1825     [ molecules ]
1826     ;molecule name   nr.
1827     Urea             1
1828     SOL              1000
1830 Ifdef statements
1831 ~~~~~~~~~~~~~~~~
1833 A very powerful feature in |Gromacs| is the use of ``#ifdef``
1834 statements in your :ref:`top` file. By making use of this
1835 statement, and associated ``#define`` statements like were
1836 seen in ``amber99.ff/forcefield.itp`` earlier,
1837 different parameters for one molecule can be used in the same
1838 :ref:`top` file. An example is given for TFE, where there is
1839 an option to use different charges on the atoms: charges derived by De
1840 Loof et al. :ref:`132 <refLoof92>` or by Van Buuren and
1841 Berendsen \ :ref:`133 <refBuuren93a>`. In fact, you can use much of the
1842 functionality of the C preprocessor, ``cpp``, because
1843 :ref:`grompp <gmx grompp>` contains similar pre-processing functions to scan
1844 the file. The way to make use of the ``#ifdef`` option is as
1845 follows:
1847 -  either use the option ``define = -DDeLoof`` in the
1848    :ref:`mdp` file (containing :ref:`grompp <gmx grompp>` input
1849    parameters), or use the line ``#define DeLoof`` early in
1850    your :ref:`top` or :ref:`itp` file; and
1852 -  put the ``#ifdef`` statements in your
1853    :ref:`top`, as shown below:
1858     ...
1862     [ atoms ]
1863     ; nr     type     resnr    residu     atom      cgnr      charge        mass
1864     #ifdef DeLoof
1865     ; Use Charges from DeLoof
1866        1        C        1        TFE        C         1        0.74        
1867        2        F        1        TFE        F         1       -0.25        
1868        3        F        1        TFE        F         1       -0.25        
1869        4        F        1        TFE        F         1       -0.25        
1870        5      CH2        1        TFE      CH2         1        0.25        
1871        6       OA        1        TFE       OA         1       -0.65        
1872        7       HO        1        TFE       HO         1        0.41        
1873     #else
1874     ; Use Charges from VanBuuren
1875        1        C        1        TFE        C         1        0.59        
1876        2        F        1        TFE        F         1       -0.2         
1877        3        F        1        TFE        F         1       -0.2         
1878        4        F        1        TFE        F         1       -0.2         
1879        5      CH2        1        TFE      CH2         1        0.26        
1880        6       OA        1        TFE       OA         1       -0.55        
1881        7       HO        1        TFE       HO         1        0.3         
1882     #endif
1884     [ bonds ]
1885     ;  ai    aj funct           c0           c1
1886         6     7     1 1.000000e-01 3.138000e+05 
1887         1     2     1 1.360000e-01 4.184000e+05 
1888         1     3     1 1.360000e-01 4.184000e+05 
1889         1     4     1 1.360000e-01 4.184000e+05 
1890         1     5     1 1.530000e-01 3.347000e+05 
1891         5     6     1 1.430000e-01 3.347000e+05 
1892     ...
1894 This mechanism is used by :ref:`pdb2gmx <gmx pdb2gmx>` to implement optional position
1895 restraints (:ref:`positionrestraint`) by ``#include``-ing an :ref:`itp` file
1896 whose contents will be meaningful only if a particular ``#define`` is set
1897 (and spelled correctly!)
1899 Topologies for free energy calculations
1900 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
1902 Free energy differences between two systems, A and B, can be calculated
1903 as described in sec. :ref:`fecalc`. Systems A and B are described by
1904 topologies consisting of the same number of molecules with the same
1905 number of atoms. Masses and non-bonded interactions can be perturbed by
1906 adding B parameters under the ``[ atoms ]`` directive. Bonded interactions can be
1907 perturbed by adding B parameters to the bonded types or the bonded
1908 interactions. The parameters that can be perturbed are listed in
1909 :numref:`Tables %s <tab-topfile1>` and :numref:`%s <tab-topfile2>`.
1910 The :math:`\lambda`-dependence of the
1911 interactions is described in section sec. :ref:`feia`. The bonded
1912 parameters that are used (on the line of the bonded interaction
1913 definition, or the ones looked up on atom types in the bonded type
1914 lists) is explained in :numref:`Table %s <tab-topfe>`. In most cases, things should
1915 work intuitively. When the A and B atom types in a bonded interaction
1916 are not all identical and parameters are not present for the B-state,
1917 either on the line or in the bonded types, :ref:`grompp <gmx grompp>` uses the A-state
1918 parameters and issues a warning. For free energy calculations, all or no
1919 parameters for topology B (:math:`\lambda = 1`) should be added on the
1920 same line, after the normal parameters, in the same order as the normal
1921 parameters. From |Gromacs| 4.6 onward, if :math:`\lambda` is treated as a
1922 vector, then the ``bonded-lambdas`` component controls all bonded terms that
1923 are not explicitly labeled as restraints. Restrain terms are controlled
1924 by the ``restraint-lambdas`` component.
1926 .. |NOT| replace:: :math:`-`
1928 .. _tab-topfe:
1930 .. table:: The bonded parameters that are used for free energy topologies,
1931            on the line of the bonded interaction definition or looked up
1932            in the bond types section based on atom types. A and B indicate the
1933            parameters used for state A and B respectively, + and |NOT| indicate
1934            the (non-)presence of parameters in the topology, x indicates that
1935            the presence has no influence.
1937            +--------------------+---------------+---------------------------------+---------+
1938            | B-state atom types | parameters    | parameters in bonded types      |         |
1939            +                    +               +-----------------+---------------+         +
1940            | all identical to   | on line       | A atom types    | B atom types  | message |
1941            +                    +-------+-------+-------+---------+-------+-------+         +
1942            | A-state atom types | A     | B     | A     | B       | A     | B     |         |
1943            +====================+=======+=======+=======+=========+=======+=======+=========+
1944            |                    | +AB   | |NOT| | x     | x       |       |       |         |
1945            |                    | +A    | +B    | x     | x       |       |       |         |
1946            | yes                | |NOT| | |NOT| | |NOT| | |NOT|   |       |       | error   |
1947            |                    | |NOT| | |NOT| | +AB   | |NOT|   |       |       |         |
1948            |                    | |NOT| | |NOT| | +A    | +B      |       |       |         |
1949            +--------------------+-------+-------+-------+---------+-------+-------+---------+
1950            |                    | +AB   | |NOT| | x     | x       | x     | x     | warning |
1951            |                    | +A    | +B    | x     | x       | x     | x     |         |
1952            |                    | |NOT| | |NOT| | |NOT| | |NOT|   | x     | x     | error   |
1953            | no                 | |NOT| | |NOT| | +AB   | |NOT|   | |NOT| | |NOT| | warning |
1954            |                    | |NOT| | |NOT| | +A    | +B      | |NOT| | |NOT| | warning |
1955            |                    | |NOT| | |NOT| | +A    | x       | +B    | |NOT| |         |
1956            |                    | |NOT| | |NOT| | +A    | x       | +     | +B    |         |
1957            +--------------------+-------+-------+-------+---------+-------+-------+---------+
1961 Below is an example of a topology which changes from 200 propanols to
1962 200 pentanes using the GROMOS-96 force field.
1966      
1967     ; Include force field parameters
1968     #include "gromos43a1.ff/forcefield.itp"
1970     [ moleculetype ]
1971     ; Name            nrexcl
1972     PropPent          3
1974     [ atoms ]
1975     ; nr type resnr residue atom cgnr  charge    mass  typeB chargeB  massB
1976       1    H    1     PROP    PH    1   0.398    1.008  CH3     0.0  15.035
1977       2   OA    1     PROP    PO    1  -0.548  15.9994  CH2     0.0  14.027
1978       3  CH2    1     PROP   PC1    1   0.150   14.027  CH2     0.0  14.027
1979       4  CH2    1     PROP   PC2    2   0.000   14.027
1980       5  CH3    1     PROP   PC3    2   0.000   15.035
1982     [ bonds ]
1983     ;  ai    aj funct    par_A  par_B 
1984         1     2     2    gb_1   gb_26
1985         2     3     2    gb_17  gb_26
1986         3     4     2    gb_26  gb_26
1987         4     5     2    gb_26
1989     [ pairs ]
1990     ;  ai    aj funct
1991         1     4     1
1992         2     5     1
1994     [ angles ]
1995     ;  ai    aj    ak funct    par_A   par_B
1996         1     2     3     2    ga_11   ga_14
1997         2     3     4     2    ga_14   ga_14
1998         3     4     5     2    ga_14   ga_14
2000     [ dihedrals ]
2001     ;  ai    aj    ak    al funct    par_A   par_B
2002         1     2     3     4     1    gd_12   gd_17
2003         2     3     4     5     1    gd_17   gd_17
2005     [ system ]
2006     ; Name
2007     Propanol to Pentane
2009     [ molecules ]
2010     ; Compound        #mols
2011     PropPent          200
2013 Atoms that are not perturbed, ``PC2`` and
2014 ``PC3``, do not need B-state parameter specifications, since
2015 the B parameters will be copied from the A parameters. Bonded
2016 interactions between atoms that are not perturbed do not need B
2017 parameter specifications, as is the case for the last bond in the
2018 example topology. Topologies using the OPLS/AA force field need no
2019 bonded parameters at all, since both the A and B parameters are
2020 determined by the atom types. Non-bonded interactions involving one or
2021 two perturbed atoms use the free-energy perturbation functional forms.
2022 Non-bonded interactions between two non-perturbed atoms use the normal
2023 functional forms. This means that when, for instance, only the charge of
2024 a particle is perturbed, its Lennard-Jones interactions will also be
2025 affected when lambda is not equal to zero or one.
2027 **Note** that this topology uses the GROMOS-96 force field, in which the
2028 bonded interactions are not determined by the atom types. The bonded
2029 interaction strings are converted by the C-preprocessor. The force-field
2030 parameter files contain lines like:
2034     #define gb_26       0.1530  7.1500e+06
2036     #define gd_17     0.000       5.86          3
2038 .. _constraintforce:
2040 Constraint forces
2041 ~~~~~~~~~~~~~~~~~
2043 | The constraint force between two atoms in one molecule can be
2044   calculated with the free energy perturbation code by adding a
2045   constraint between the two atoms, with a different length in the A and
2046   B topology. When the B length is 1 nm longer than the A length and
2047   lambda is kept constant at zero, the derivative of the Hamiltonian
2048   with respect to lambda is the constraint force. For constraints
2049   between molecules, the pull code can be used, see sec. :ref:`pull`.
2050   Below is an example for calculating the constraint force at 0.7 nm
2051   between two methanes in water, by combining the two methanes into one
2052   “molecule.” **Note** that the definition of a “molecule” in |Gromacs|
2053   does not necessarily correspond to the chemical definition of a
2054   molecule. In |Gromacs|, a “molecule” can be defined as any group of
2055   atoms that one wishes to consider simultaneously. The added constraint
2056   is of function type 2, which means that it is not used for generating
2057   exclusions (see sec. :ref:`excl`). Note that the constraint free energy
2058   term is included in the derivative term, and is specifically included
2059   in the ``bonded-lambdas`` component. However, the free energy for changing
2060   constraints is *not* included in the potential energy differences used
2061   for BAR and MBAR, as this requires reevaluating the energy at each of
2062   the constraint components. This functionality is planned for later
2063   versions.
2067     ; Include force-field parameters
2068     #include "gromos43a1.ff/forcefield.itp"
2070     [ moleculetype ]
2071     ; Name            nrexcl
2072     Methanes               1
2074     [ atoms ]
2075     ; nr   type   resnr  residu   atom    cgnr     charge    mass
2076        1    CH4     1     CH4      C1       1          0    16.043
2077        2    CH4     1     CH4      C2       2          0    16.043
2078     [ constraints ]
2079     ;  ai    aj funct   length_A  length_B
2080         1     2     2        0.7       1.7
2082     #include "gromos43a1.ff/spc.itp"
2084     [ system ]
2085     ; Name
2086     Methanes in Water
2088     [ molecules ]
2089     ; Compound        #mols
2090     Methanes              1
2091     SOL                2002
2093 Coordinate file
2094 ~~~~~~~~~~~~~~~
2096 Files with the :ref:`gro` file extension contain a molecular
2097 structure in GROMOS-87 format. A sample piece is included below:
2101     MD of 2 waters, reformat step, PA aug-91
2102         6
2103         1WATER  OW1    1   0.126   1.624   1.679  0.1227 -0.0580  0.0434
2104         1WATER  HW2    2   0.190   1.661   1.747  0.8085  0.3191 -0.7791
2105         1WATER  HW3    3   0.177   1.568   1.613 -0.9045 -2.6469  1.3180
2106         2WATER  OW1    4   1.275   0.053   0.622  0.2519  0.3140 -0.1734
2107         2WATER  HW2    5   1.337   0.002   0.680 -1.0641 -1.1349  0.0257
2108         2WATER  HW3    6   1.326   0.120   0.568  1.9427 -0.8216 -0.0244
2109        1.82060   1.82060   1.82060
2111 This format is fixed, *i.e.* all columns are in a fixed position. If you
2112 want to read such a file in your own program without using the |Gromacs|
2113 libraries you can use the following formats:
2115 **C-format:**
2116 ``“%5i%5s%5s%5i%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f”``
2118 Or to be more precise, with title *etc.* it looks like this:
2122       "%s\n", Title
2123       "%5d\n", natoms
2124       for (i=0; (i<natoms); i++) {
2125         "%5d%-5s%5s%5d%8.3f%8.3f%8.3f%8.4f%8.4f%8.4f\n",
2126           residuenr,residuename,atomname,atomnr,x,y,z,vx,vy,vz
2127       }
2128       "%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f%10.5f\n",
2129         box[X][X],box[Y][Y],box[Z][Z],
2130         box[X][Y],box[X][Z],box[Y][X],box[Y][Z],box[Z][X],box[Z][Y]
2132 **Fortran format:**
2133 ``(i5,2a5,i5,3f8.3,3f8.4)``
2135 So ``confin.gro`` is the |Gromacs| coordinate file and is
2136 almost the same as the GROMOS-87 file (for GROMOS users: when used with
2137 ``ntx=7``). The only difference is the box for which |Gromacs|
2138 uses a tensor, not a vector.
2140 .. _fforganization:
2142 Force field organization
2143 ------------------------
2145 .. _fffiles:
2147 Force-field files
2148 ~~~~~~~~~~~~~~~~~
2150 Many force fields are available by default. Force fields are detected by
2151 the presence of ``<name>.ff`` directories in the
2152 ``$GMXLIB/share/gromacs/top`` sub-directory and/or the
2153 working directory. The information regarding the location of the force
2154 field files is printed by :ref:`pdb2gmx <gmx pdb2gmx>` so you can easily keep
2155 track of which version of a force field is being called, in case you
2156 have made modifications in one location or another. The force fields
2157 included with |Gromacs| are:
2159 -  AMBER03 protein, nucleic AMBER94 (Duan et al., J. Comp. Chem. 24,
2160    1999-2012, 2003)
2162 -  AMBER94 force field (Cornell et al., JACS 117, 5179-5197, 1995)
2164 -  AMBER96 protein, nucleic AMBER94 (Kollman et al., Acc. Chem. Res. 29,
2165    461-469, 1996)
2167 -  AMBER99 protein, nucleic AMBER94 (Wang et al., J. Comp. Chem. 21,
2168    1049-1074, 2000)
2170 -  AMBER99SB protein, nucleic AMBER94 (Hornak et al., Proteins 65,
2171    712-725, 2006)
2173 -  AMBER99SB-ILDN protein, nucleic AMBER94 (Lindorff-Larsen et al.,
2174    Proteins 78, 1950-58, 2010)
2176 -  AMBERGS force field (Garcia & Sanbonmatsu, PNAS 99, 2782-2787, 2002)
2178 -  CHARMM27 all-atom force field (CHARM22 plus CMAP for proteins)
2180 -  GROMOS96 43a1 force field
2182 -  GROMOS96 43a2 force field (improved alkane dihedrals)
2184 -  GROMOS96 45a3 force field (Schuler JCC 2001 22 1205)
2186 -  GROMOS96 53a5 force field (JCC 2004 vol 25 pag 1656)
2188 -  GROMOS96 53a6 force field (JCC 2004 vol 25 pag 1656)
2190 -  GROMOS96 54a7 force field (Eur. Biophys. J. (2011), 40,, 843-856,
2191    DOI: 10.1007/s00249-011-0700-9)
2193 -  OPLS-AA/L all-atom force field (2001 aminoacid dihedrals)
2195 A force field is included at the beginning of a topology file with an
2196 ``#include`` statement followed by
2197 ``<name>.ff/forcefield.itp``. This statement includes the
2198 force-field file, which, in turn, may include other force-field files.
2199 All the force fields are organized in the same way. An example of the
2200 ``amber99.ff/forcefield.itp`` was shown in
2201 :ref:`topfile`.
2203 For each force field, there several files which are only used by
2204 :ref:`pdb2gmx <gmx pdb2gmx>`. These are: residue databases
2205 (:ref:`rtp`) the hydrogen
2206 database (:ref:`hdb`), two
2207 termini databases (``.n.tdb`` and ``.c.tdb``,
2208 see ) and the atom type database
2209 (:ref:`atp`), which
2210 contains only the masses. Other optional files are described in sec. :ref:`pdb2gmxfiles`.
2212 Changing force-field parameters
2213 ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
2215 If one wants to change the parameters of few bonded interactions in a
2216 molecule, this is most easily accomplished by typing the parameters
2217 behind the definition of the bonded interaction directly in the
2218 :ref:`top` file under the ``[ moleculetype ]``
2219 section (see :ref:`topfile` for the format and units).
2220 If one wants to change the parameters for all instances of a
2221 certain interaction one can change them in the force-field file or add a
2222 new ``[ ???types ]`` section after including the force
2223 field. When parameters for a certain interaction are defined multiple
2224 times, the last definition is used. As of |Gromacs| version 3.1.3, a
2225 warning is generated when parameters are redefined with a different
2226 value. Changing the Lennard-Jones parameters of an atom type is not
2227 recommended, because in the GROMOS force fields the Lennard-Jones
2228 parameters for several combinations of atom types are not generated
2229 according to the standard combination rules. Such combinations (and
2230 possibly others that do follow the combination rules) are defined in the
2231 ``[ nonbond_params ]`` section, and changing the
2232 Lennard-Jones parameters of an atom type has no effect on these
2233 combinations.
2235 Adding atom types
2236 ~~~~~~~~~~~~~~~~~
2238 As of |Gromacs| version 3.1.3, atom types can be added in an extra
2239 ``[ atomtypes ]`` section after the inclusion of the
2240 normal force field. After the definition of the new atom type(s),
2241 additional non-bonded and pair parameters can be defined. In pre-3.1.3
2242 versions of |Gromacs|, the new atom types needed to be added in the
2243 ``[ atomtypes ]`` section of the force-field files, because
2244 all non-bonded parameters above the last ``[ atomtypes ]``
2245 section would be overwritten using the standard combination rules.