2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
37 /* This file is completely threadsafe - keep it that way! */
46 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
47 #include "gromacs/gmxpreprocess/topio.h"
48 #include "gromacs/gmxpreprocess/toputil.h"
49 #include "gromacs/legacyheaders/names.h"
50 #include "gromacs/legacyheaders/typedefs.h"
51 #include "gromacs/legacyheaders/types/ifunc.h"
52 #include "gromacs/math/units.h"
53 #include "gromacs/math/utilities.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/utility/fatalerror.h"
56 #include "gromacs/utility/smalloc.h"
58 static int round_check(real r
, int limit
, int ftype
, const char *name
)
71 if (r
-i
> 0.01 || r
-i
< -0.01)
73 gmx_fatal(FARGS
, "A non-integer value (%f) was supplied for '%s' in %s",
74 r
, name
, interaction_function
[ftype
].longname
);
79 gmx_fatal(FARGS
, "Value of '%s' in %s is %d, which is smaller than the minimum of %d",
80 name
, interaction_function
[ftype
].longname
, i
, limit
);
86 static void set_ljparams(int comb
, double reppow
, double v
, double w
,
89 if (comb
== eCOMB_ARITHMETIC
|| comb
== eCOMB_GEOM_SIG_EPS
)
93 *c6
= 4*w
*std::pow(v
, 6.0);
94 *c12
= 4*w
*std::pow(v
, reppow
);
98 /* Interpret negative sigma as c6=0 and c12 with -sigma */
100 *c12
= 4*w
*std::pow(-v
, reppow
);
110 /* A return value of 0 means parameters were assigned successfully,
111 * returning -1 means this is an all-zero interaction that should not be added.
114 assign_param(t_functype ftype
, t_iparams
*newparam
,
115 real old
[MAXFORCEPARAM
], int comb
, double reppow
)
118 gmx_bool all_param_zero
= TRUE
;
121 for (j
= 0; (j
< MAXFORCEPARAM
); j
++)
123 newparam
->generic
.buf
[j
] = 0.0;
124 /* If all parameters are zero we might not add some interaction types (selected below).
125 * We cannot apply this to ALL interactions, since many have valid reasons for having
126 * zero parameters (e.g. an index to a Cmap interaction, or LJ parameters), but
127 * we use it for angles and torsions that are typically generated automatically.
129 all_param_zero
= (all_param_zero
== TRUE
) && fabs(old
[j
]) < GMX_REAL_MIN
;
132 if (all_param_zero
== TRUE
)
134 if (IS_ANGLE(ftype
) || IS_RESTRAINT_TYPE(ftype
) || ftype
== F_IDIHS
||
135 ftype
== F_PDIHS
|| ftype
== F_PIDIHS
|| ftype
== F_RBDIHS
|| ftype
== F_FOURDIHS
)
144 /* Post processing of input data: store cosine iso angle itself */
145 newparam
->harmonic
.rA
= cos(old
[0]*DEG2RAD
);
146 newparam
->harmonic
.krA
= old
[1];
147 newparam
->harmonic
.rB
= cos(old
[2]*DEG2RAD
);
148 newparam
->harmonic
.krB
= old
[3];
151 /* Post processing of input data: store square of length itself */
152 newparam
->harmonic
.rA
= sqr(old
[0]);
153 newparam
->harmonic
.krA
= old
[1];
154 newparam
->harmonic
.rB
= sqr(old
[2]);
155 newparam
->harmonic
.krB
= old
[3];
158 newparam
->fene
.bm
= old
[0];
159 newparam
->fene
.kb
= old
[1];
162 newparam
->restraint
.lowA
= old
[0];
163 newparam
->restraint
.up1A
= old
[1];
164 newparam
->restraint
.up2A
= old
[2];
165 newparam
->restraint
.kA
= old
[3];
166 newparam
->restraint
.lowB
= old
[4];
167 newparam
->restraint
.up1B
= old
[5];
168 newparam
->restraint
.up2B
= old
[6];
169 newparam
->restraint
.kB
= old
[7];
175 newparam
->tab
.table
= round_check(old
[0], 0, ftype
, "table index");
176 newparam
->tab
.kA
= old
[1];
177 newparam
->tab
.kB
= old
[3];
179 case F_CROSS_BOND_BONDS
:
180 newparam
->cross_bb
.r1e
= old
[0];
181 newparam
->cross_bb
.r2e
= old
[1];
182 newparam
->cross_bb
.krr
= old
[2];
184 case F_CROSS_BOND_ANGLES
:
185 newparam
->cross_ba
.r1e
= old
[0];
186 newparam
->cross_ba
.r2e
= old
[1];
187 newparam
->cross_ba
.r3e
= old
[2];
188 newparam
->cross_ba
.krt
= old
[3];
191 newparam
->u_b
.thetaA
= old
[0];
192 newparam
->u_b
.kthetaA
= old
[1];
193 newparam
->u_b
.r13A
= old
[2];
194 newparam
->u_b
.kUBA
= old
[3];
195 newparam
->u_b
.thetaB
= old
[4];
196 newparam
->u_b
.kthetaB
= old
[5];
197 newparam
->u_b
.r13B
= old
[6];
198 newparam
->u_b
.kUBB
= old
[7];
200 case F_QUARTIC_ANGLES
:
201 newparam
->qangle
.theta
= old
[0];
202 for (i
= 0; i
< 5; i
++)
204 newparam
->qangle
.c
[i
] = old
[i
+1];
207 case F_LINEAR_ANGLES
:
208 newparam
->linangle
.aA
= old
[0];
209 newparam
->linangle
.klinA
= old
[1];
210 newparam
->linangle
.aB
= old
[2];
211 newparam
->linangle
.klinB
= old
[3];
217 newparam
->harmonic
.rA
= old
[0];
218 newparam
->harmonic
.krA
= old
[1];
219 newparam
->harmonic
.rB
= old
[2];
220 newparam
->harmonic
.krB
= old
[3];
223 newparam
->harmonic
.rA
= old
[0];
224 newparam
->harmonic
.krA
= old
[1];
227 newparam
->morse
.b0A
= old
[0];
228 newparam
->morse
.cbA
= old
[1];
229 newparam
->morse
.betaA
= old
[2];
230 newparam
->morse
.b0B
= old
[3];
231 newparam
->morse
.cbB
= old
[4];
232 newparam
->morse
.betaB
= old
[5];
235 newparam
->cubic
.b0
= old
[0];
236 newparam
->cubic
.kb
= old
[1];
237 newparam
->cubic
.kcub
= old
[2];
242 newparam
->polarize
.alpha
= old
[0];
245 newparam
->anharm_polarize
.alpha
= old
[0];
246 newparam
->anharm_polarize
.drcut
= old
[1];
247 newparam
->anharm_polarize
.khyp
= old
[2];
250 newparam
->wpol
.al_x
= old
[0];
251 newparam
->wpol
.al_y
= old
[1];
252 newparam
->wpol
.al_z
= old
[2];
253 newparam
->wpol
.rOH
= old
[3];
254 newparam
->wpol
.rHH
= old
[4];
255 newparam
->wpol
.rOD
= old
[5];
258 newparam
->thole
.a
= old
[0];
259 newparam
->thole
.alpha1
= old
[1];
260 newparam
->thole
.alpha2
= old
[2];
261 if ((old
[1] > 0) && (old
[2] > 0))
263 newparam
->thole
.rfac
= old
[0]*std::pow(old
[1]*old
[2], static_cast<real
>(-1.0/6.0));
267 newparam
->thole
.rfac
= 1;
271 newparam
->bham
.a
= old
[0];
272 newparam
->bham
.b
= old
[1];
273 newparam
->bham
.c
= old
[2];
276 set_ljparams(comb
, reppow
, old
[0], old
[1], &newparam
->lj14
.c6A
, &newparam
->lj14
.c12A
);
277 set_ljparams(comb
, reppow
, old
[2], old
[3], &newparam
->lj14
.c6B
, &newparam
->lj14
.c12B
);
280 newparam
->ljc14
.fqq
= old
[0];
281 newparam
->ljc14
.qi
= old
[1];
282 newparam
->ljc14
.qj
= old
[2];
283 set_ljparams(comb
, reppow
, old
[3], old
[4], &newparam
->ljc14
.c6
, &newparam
->ljc14
.c12
);
286 newparam
->ljcnb
.qi
= old
[0];
287 newparam
->ljcnb
.qj
= old
[1];
288 set_ljparams(comb
, reppow
, old
[2], old
[3], &newparam
->ljcnb
.c6
, &newparam
->ljcnb
.c12
);
291 set_ljparams(comb
, reppow
, old
[0], old
[1], &newparam
->lj
.c6
, &newparam
->lj
.c12
);
297 newparam
->pdihs
.phiA
= old
[0];
298 newparam
->pdihs
.cpA
= old
[1];
300 /* Change 20100720: Amber occasionally uses negative multiplicities (mathematically OK),
301 * so I have changed the lower limit to -99 /EL
303 newparam
->pdihs
.phiB
= old
[3];
304 newparam
->pdihs
.cpB
= old
[4];
305 /* If both force constants are zero there is no interaction. Return -1 to signal
306 * this entry should NOT be added.
308 if (fabs(newparam
->pdihs
.cpA
) < GMX_REAL_MIN
&& fabs(newparam
->pdihs
.cpB
) < GMX_REAL_MIN
)
313 newparam
->pdihs
.mult
= round_check(old
[2], -99, ftype
, "multiplicity");
317 newparam
->pdihs
.phiA
= old
[0];
318 newparam
->pdihs
.cpA
= old
[1];
321 newparam
->posres
.fcA
[XX
] = old
[0];
322 newparam
->posres
.fcA
[YY
] = old
[1];
323 newparam
->posres
.fcA
[ZZ
] = old
[2];
324 newparam
->posres
.fcB
[XX
] = old
[3];
325 newparam
->posres
.fcB
[YY
] = old
[4];
326 newparam
->posres
.fcB
[ZZ
] = old
[5];
327 newparam
->posres
.pos0A
[XX
] = old
[6];
328 newparam
->posres
.pos0A
[YY
] = old
[7];
329 newparam
->posres
.pos0A
[ZZ
] = old
[8];
330 newparam
->posres
.pos0B
[XX
] = old
[9];
331 newparam
->posres
.pos0B
[YY
] = old
[10];
332 newparam
->posres
.pos0B
[ZZ
] = old
[11];
335 newparam
->fbposres
.geom
= round_check(old
[0], 0, ftype
, "geometry");
336 if (!(newparam
->fbposres
.geom
> efbposresZERO
&& newparam
->fbposres
.geom
< efbposresNR
))
338 gmx_fatal(FARGS
, "Invalid geometry for flat-bottomed position restraint.\n"
339 "Expected number between 1 and %d. Found %d\n", efbposresNR
-1,
340 newparam
->fbposres
.geom
);
342 newparam
->fbposres
.r
= old
[1];
343 newparam
->fbposres
.k
= old
[2];
344 newparam
->fbposres
.pos0
[XX
] = old
[3];
345 newparam
->fbposres
.pos0
[YY
] = old
[4];
346 newparam
->fbposres
.pos0
[ZZ
] = old
[5];
349 newparam
->disres
.label
= round_check(old
[0], 0, ftype
, "label");
350 newparam
->disres
.type
= round_check(old
[1], 1, ftype
, "type'");
351 newparam
->disres
.low
= old
[2];
352 newparam
->disres
.up1
= old
[3];
353 newparam
->disres
.up2
= old
[4];
354 newparam
->disres
.kfac
= old
[5];
357 newparam
->orires
.ex
= round_check(old
[0], 1, ftype
, "experiment") - 1;
358 newparam
->orires
.label
= round_check(old
[1], 1, ftype
, "label");
359 newparam
->orires
.power
= round_check(old
[2], 0, ftype
, "power");
360 newparam
->orires
.c
= old
[3];
361 newparam
->orires
.obs
= old
[4];
362 newparam
->orires
.kfac
= old
[5];
365 newparam
->dihres
.phiA
= old
[0];
366 newparam
->dihres
.dphiA
= old
[1];
367 newparam
->dihres
.kfacA
= old
[2];
368 newparam
->dihres
.phiB
= old
[3];
369 newparam
->dihres
.dphiB
= old
[4];
370 newparam
->dihres
.kfacB
= old
[5];
373 for (i
= 0; (i
< NR_RBDIHS
); i
++)
375 newparam
->rbdihs
.rbcA
[i
] = old
[i
];
376 newparam
->rbdihs
.rbcB
[i
] = old
[NR_RBDIHS
+i
];
380 for (i
= 0; (i
< NR_CBTDIHS
); i
++)
382 newparam
->cbtdihs
.cbtcA
[i
] = old
[i
];
386 /* Read the dihedral parameters to temporary arrays,
387 * and convert them to the computationally faster
388 * Ryckaert-Bellemans form.
390 /* Use conversion formula for OPLS to Ryckaert-Bellemans: */
391 newparam
->rbdihs
.rbcA
[0] = old
[1]+0.5*(old
[0]+old
[2]);
392 newparam
->rbdihs
.rbcA
[1] = 0.5*(3.0*old
[2]-old
[0]);
393 newparam
->rbdihs
.rbcA
[2] = 4.0*old
[3]-old
[1];
394 newparam
->rbdihs
.rbcA
[3] = -2.0*old
[2];
395 newparam
->rbdihs
.rbcA
[4] = -4.0*old
[3];
396 newparam
->rbdihs
.rbcA
[5] = 0.0;
398 newparam
->rbdihs
.rbcB
[0] = old
[NR_FOURDIHS
+1]+0.5*(old
[NR_FOURDIHS
+0]+old
[NR_FOURDIHS
+2]);
399 newparam
->rbdihs
.rbcB
[1] = 0.5*(3.0*old
[NR_FOURDIHS
+2]-old
[NR_FOURDIHS
+0]);
400 newparam
->rbdihs
.rbcB
[2] = 4.0*old
[NR_FOURDIHS
+3]-old
[NR_FOURDIHS
+1];
401 newparam
->rbdihs
.rbcB
[3] = -2.0*old
[NR_FOURDIHS
+2];
402 newparam
->rbdihs
.rbcB
[4] = -4.0*old
[NR_FOURDIHS
+3];
403 newparam
->rbdihs
.rbcB
[5] = 0.0;
407 newparam
->constr
.dA
= old
[0];
408 newparam
->constr
.dB
= old
[1];
411 newparam
->settle
.doh
= old
[0];
412 newparam
->settle
.dhh
= old
[1];
420 newparam
->vsite
.a
= old
[0];
421 newparam
->vsite
.b
= old
[1];
422 newparam
->vsite
.c
= old
[2];
423 newparam
->vsite
.d
= old
[3];
424 newparam
->vsite
.e
= old
[4];
425 newparam
->vsite
.f
= old
[5];
428 newparam
->vsite
.a
= old
[1] * cos(DEG2RAD
* old
[0]);
429 newparam
->vsite
.b
= old
[1] * sin(DEG2RAD
* old
[0]);
430 newparam
->vsite
.c
= old
[2];
431 newparam
->vsite
.d
= old
[3];
432 newparam
->vsite
.e
= old
[4];
433 newparam
->vsite
.f
= old
[5];
436 newparam
->vsiten
.n
= round_check(old
[0], 1, ftype
, "number of atoms");
437 newparam
->vsiten
.a
= old
[1];
440 newparam
->cmap
.cmapA
= static_cast<int>(old
[0]);
441 newparam
->cmap
.cmapB
= static_cast<int>(old
[1]);
446 newparam
->gb
.sar
= old
[0];
447 newparam
->gb
.st
= old
[1];
448 newparam
->gb
.pi
= old
[2];
449 newparam
->gb
.gbr
= old
[3];
450 newparam
->gb
.bmlt
= old
[4];
453 gmx_fatal(FARGS
, "unknown function type %d in %s line %d",
454 ftype
, __FILE__
, __LINE__
);
459 static int enter_params(gmx_ffparams_t
*ffparams
, t_functype ftype
,
460 real forceparams
[MAXFORCEPARAM
], int comb
, real reppow
,
461 int start
, gmx_bool bAppend
)
467 if ( (rc
= assign_param(ftype
, &newparam
, forceparams
, comb
, reppow
)) < 0)
469 /* -1 means this interaction is all-zero and should not be added */
475 for (type
= start
; (type
< ffparams
->ntypes
); type
++)
477 if (ffparams
->functype
[type
] == ftype
)
481 /* Occasionally, the way the 1-3 reference distance is
482 * computed can lead to non-binary-identical results, but I
484 if ((gmx_within_tol(newparam
.gb
.sar
, ffparams
->iparams
[type
].gb
.sar
, 1e-6)) &&
485 (gmx_within_tol(newparam
.gb
.st
, ffparams
->iparams
[type
].gb
.st
, 1e-6)) &&
486 (gmx_within_tol(newparam
.gb
.pi
, ffparams
->iparams
[type
].gb
.pi
, 1e-6)) &&
487 (gmx_within_tol(newparam
.gb
.gbr
, ffparams
->iparams
[type
].gb
.gbr
, 1e-6)) &&
488 (gmx_within_tol(newparam
.gb
.bmlt
, ffparams
->iparams
[type
].gb
.bmlt
, 1e-6)))
495 if (memcmp(&newparam
, &ffparams
->iparams
[type
], (size_t)sizeof(newparam
)) == 0)
505 type
= ffparams
->ntypes
;
509 fprintf(debug
, "copying newparam to ffparams->iparams[%d] (ntypes=%d)\n",
510 type
, ffparams
->ntypes
);
512 memcpy(&ffparams
->iparams
[type
], &newparam
, (size_t)sizeof(newparam
));
515 ffparams
->functype
[type
] = ftype
;
520 static void append_interaction(t_ilist
*ilist
,
521 int type
, int nral
, atom_id a
[MAXATOMLIST
])
528 ilist
->iatoms
[where1
++] = type
;
529 for (i
= 0; (i
< nral
); i
++)
531 ilist
->iatoms
[where1
++] = a
[i
];
535 static void enter_function(t_params
*p
, t_functype ftype
, int comb
, real reppow
,
536 gmx_ffparams_t
*ffparams
, t_ilist
*il
,
538 gmx_bool bNB
, gmx_bool bAppend
)
540 int k
, type
, nr
, nral
, delta
, start
;
542 start
= ffparams
->ntypes
;
545 for (k
= 0; k
< nr
; k
++)
547 if (*maxtypes
<= ffparams
->ntypes
)
550 srenew(ffparams
->functype
, *maxtypes
);
551 srenew(ffparams
->iparams
, *maxtypes
);
554 fprintf(debug
, "%s, line %d: srenewed idef->functype and idef->iparams to %d\n",
555 __FILE__
, __LINE__
, *maxtypes
);
558 type
= enter_params(ffparams
, ftype
, p
->param
[k
].c
, comb
, reppow
, start
, bAppend
);
559 /* Type==-1 is used as a signal that this interaction is all-zero and should not be added. */
560 if (!bNB
&& type
>= 0)
564 srenew(il
->iatoms
, il
->nr
+delta
);
565 append_interaction(il
, type
, nral
, p
->param
[k
].a
);
570 void convert_params(int atnr
, t_params nbtypes
[],
571 t_molinfo
*mi
, t_molinfo
*intermolecular_interactions
,
572 int comb
, double reppow
, real fudgeQQ
,
583 ffp
= &mtop
->ffparams
;
586 ffp
->functype
= NULL
;
588 ffp
->reppow
= reppow
;
590 enter_function(&(nbtypes
[F_LJ
]), (t_functype
)F_LJ
, comb
, reppow
, ffp
, NULL
,
591 &maxtypes
, TRUE
, TRUE
);
592 enter_function(&(nbtypes
[F_BHAM
]), (t_functype
)F_BHAM
, comb
, reppow
, ffp
, NULL
,
593 &maxtypes
, TRUE
, TRUE
);
595 for (mt
= 0; mt
< mtop
->nmoltype
; mt
++)
597 molt
= &mtop
->moltype
[mt
];
598 for (i
= 0; (i
< F_NRE
); i
++)
600 molt
->ilist
[i
].nr
= 0;
601 molt
->ilist
[i
].iatoms
= NULL
;
603 plist
= mi
[mt
].plist
;
605 flags
= interaction_function
[i
].flags
;
606 if ((i
!= F_LJ
) && (i
!= F_BHAM
) && ((flags
& IF_BOND
) ||
607 (flags
& IF_VSITE
) ||
608 (flags
& IF_CONSTRAINT
)))
610 enter_function(&(plist
[i
]), (t_functype
)i
, comb
, reppow
,
611 ffp
, &molt
->ilist
[i
],
612 &maxtypes
, FALSE
, (i
== F_POSRES
|| i
== F_FBPOSRES
));
617 mtop
->bIntermolecularInteractions
= FALSE
;
618 if (intermolecular_interactions
!= NULL
)
620 /* Process the intermolecular interaction list */
621 snew(mtop
->intermolecular_ilist
, F_NRE
);
623 for (i
= 0; (i
< F_NRE
); i
++)
625 mtop
->intermolecular_ilist
[i
].nr
= 0;
626 mtop
->intermolecular_ilist
[i
].iatoms
= NULL
;
628 plist
= intermolecular_interactions
->plist
;
632 flags
= interaction_function
[i
].flags
;
633 /* For intermolecular interactions we (currently)
634 * only support potentials.
635 * Constraints and virtual sites would be possible,
636 * but require a lot of extra (bug-prone) code.
638 if (!(flags
& IF_BOND
))
640 gmx_fatal(FARGS
, "The intermolecular_interaction section may only contain bonded potentials");
642 else if (NRAL(i
) == 1) /* e.g. position restraints */
644 gmx_fatal(FARGS
, "Single atom interactions don't make sense in the intermolecular_interaction section, you can put them in the moleculetype section");
646 else if (flags
& IF_CHEMBOND
)
648 gmx_fatal(FARGS
, "The intermolecular_interaction can not contain chemically bonding interactions");
652 enter_function(&(plist
[i
]), (t_functype
)i
, comb
, reppow
,
653 ffp
, &mtop
->intermolecular_ilist
[i
],
654 &maxtypes
, FALSE
, FALSE
);
656 mtop
->bIntermolecularInteractions
= TRUE
;
661 if (!mtop
->bIntermolecularInteractions
)
663 sfree(mtop
->intermolecular_ilist
);
669 fprintf(debug
, "%s, line %d: There are %d functypes in idef\n",
670 __FILE__
, __LINE__
, ffp
->ntypes
);
673 ffp
->fudgeQQ
= fudgeQQ
;