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 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
38 /* This file is completely threadsafe - keep it that way! */
49 #include "gmx_fatal.h"
54 #include "gpp_atomtype.h"
57 static int round_check(real r
,int limit
,int ftype
,const char *name
)
66 if (r
-i
> 0.01 || r
-i
< -0.01)
67 gmx_fatal(FARGS
,"A non-integer value (%f) was supplied for '%s' in %s",
68 r
,name
,interaction_function
[ftype
].longname
);
71 gmx_fatal(FARGS
,"Value of '%s' in %s is %d, which is smaller than the minimum of %d",
72 name
,interaction_function
[ftype
].longname
,i
,limit
);
77 static void set_ljparams(int comb
,double reppow
,real v
,real w
,
80 if (comb
== eCOMB_ARITHMETIC
|| comb
== eCOMB_GEOM_SIG_EPS
) {
83 *c12
= 4*w
*pow(v
,reppow
);
85 /* Interpret negative sigma as c6=0 and c12 with -sigma */
87 *c12
= 4*w
*pow(-v
,reppow
);
95 /* A return value of 0 means parameters were assigned successfully,
96 * returning -1 means this is an all-zero interaction that should not be added.
99 assign_param(t_functype ftype
,t_iparams
*newparam
,
100 real old
[MAXFORCEPARAM
],int comb
,double reppow
)
104 gmx_bool all_param_zero
=TRUE
;
107 for(j
=0; (j
<MAXFORCEPARAM
); j
++)
109 newparam
->generic
.buf
[j
]=0.0;
110 /* If all parameters are zero we might not add some interaction types (selected below).
111 * We cannot apply this to ALL interactions, since many have valid reasons for having
112 * zero parameters (e.g. an index to a Cmap interaction, or LJ parameters), but
113 * we use it for angles and torsions that are typically generated automatically.
115 all_param_zero
= (all_param_zero
==TRUE
) && fabs(old
[j
])<GMX_REAL_MIN
;
118 if(all_param_zero
==TRUE
)
120 if(IS_ANGLE(ftype
) || IS_RESTRAINT_TYPE(ftype
) || ftype
==F_IDIHS
||
121 ftype
==F_PDIHS
|| ftype
==F_PIDIHS
|| ftype
==F_RBDIHS
|| ftype
==F_FOURDIHS
)
129 /* Post processing of input data: store cosine iso angle itself */
130 newparam
->harmonic
.rA
=cos(old
[0]*DEG2RAD
);
131 newparam
->harmonic
.krA
=old
[1];
132 newparam
->harmonic
.rB
=cos(old
[2]*DEG2RAD
);
133 newparam
->harmonic
.krB
=old
[3];
136 /* Post processing of input data: store square of length itself */
137 newparam
->harmonic
.rA
=sqr(old
[0]);
138 newparam
->harmonic
.krA
=old
[1];
139 newparam
->harmonic
.rB
=sqr(old
[2]);
140 newparam
->harmonic
.krB
=old
[3];
143 newparam
->fene
.bm
=old
[0];
144 newparam
->fene
.kb
=old
[1];
147 newparam
->restraint
.lowA
= old
[0];
148 newparam
->restraint
.up1A
= old
[1];
149 newparam
->restraint
.up2A
= old
[2];
150 newparam
->restraint
.kA
= old
[3];
151 newparam
->restraint
.lowB
= old
[4];
152 newparam
->restraint
.up1B
= old
[5];
153 newparam
->restraint
.up2B
= old
[6];
154 newparam
->restraint
.kB
= old
[7];
160 newparam
->tab
.table
= round_check(old
[0],0,ftype
,"table index");
161 newparam
->tab
.kA
= old
[1];
162 newparam
->tab
.kB
= old
[3];
164 case F_CROSS_BOND_BONDS
:
165 newparam
->cross_bb
.r1e
=old
[0];
166 newparam
->cross_bb
.r2e
=old
[1];
167 newparam
->cross_bb
.krr
=old
[2];
169 case F_CROSS_BOND_ANGLES
:
170 newparam
->cross_ba
.r1e
=old
[0];
171 newparam
->cross_ba
.r2e
=old
[1];
172 newparam
->cross_ba
.r3e
=old
[2];
173 newparam
->cross_ba
.krt
=old
[3];
176 newparam
->u_b
.thetaA
=old
[0];
177 newparam
->u_b
.kthetaA
=old
[1];
178 newparam
->u_b
.r13A
=old
[2];
179 newparam
->u_b
.kUBA
=old
[3];
180 newparam
->u_b
.thetaB
=old
[4];
181 newparam
->u_b
.kthetaB
=old
[5];
182 newparam
->u_b
.r13B
=old
[6];
183 newparam
->u_b
.kUBB
=old
[7];
185 case F_QUARTIC_ANGLES
:
186 newparam
->qangle
.theta
=old
[0];
188 newparam
->qangle
.c
[i
]=old
[i
+1];
190 case F_LINEAR_ANGLES
:
191 newparam
->linangle
.aA
= old
[0];
192 newparam
->linangle
.klinA
= old
[1];
193 newparam
->linangle
.aB
= old
[2];
194 newparam
->linangle
.klinB
= old
[3];
200 newparam
->harmonic
.rA
=old
[0];
201 newparam
->harmonic
.krA
=old
[1];
202 newparam
->harmonic
.rB
=old
[2];
203 newparam
->harmonic
.krB
=old
[3];
206 newparam
->morse
.b0A
=old
[0];
207 newparam
->morse
.cbA
=old
[1];
208 newparam
->morse
.betaA
=old
[2];
209 newparam
->morse
.b0B
=old
[3];
210 newparam
->morse
.cbB
=old
[4];
211 newparam
->morse
.betaB
=old
[5];
214 newparam
->cubic
.b0
=old
[0];
215 newparam
->cubic
.kb
=old
[1];
216 newparam
->cubic
.kcub
=old
[2];
221 newparam
->polarize
.alpha
= old
[0];
224 newparam
->anharm_polarize
.alpha
= old
[0];
225 newparam
->anharm_polarize
.drcut
= old
[1];
226 newparam
->anharm_polarize
.khyp
= old
[2];
229 newparam
->wpol
.al_x
=old
[0];
230 newparam
->wpol
.al_y
=old
[1];
231 newparam
->wpol
.al_z
=old
[2];
232 newparam
->wpol
.rOH
=old
[3];
233 newparam
->wpol
.rHH
=old
[4];
234 newparam
->wpol
.rOD
=old
[5];
237 newparam
->thole
.a
= old
[0];
238 newparam
->thole
.alpha1
= old
[1];
239 newparam
->thole
.alpha2
= old
[2];
240 if ((old
[1] > 0) && (old
[2] > 0))
241 newparam
->thole
.rfac
= old
[0]*pow(old
[1]*old
[2],-1.0/6.0);
243 newparam
->thole
.rfac
= 1;
246 newparam
->bham
.a
= old
[0];
247 newparam
->bham
.b
= old
[1];
248 newparam
->bham
.c
= old
[2];
251 set_ljparams(comb
,reppow
,old
[0],old
[1],&newparam
->lj14
.c6A
,&newparam
->lj14
.c12A
);
252 set_ljparams(comb
,reppow
,old
[2],old
[3],&newparam
->lj14
.c6B
,&newparam
->lj14
.c12B
);
255 newparam
->ljc14
.fqq
= old
[0];
256 newparam
->ljc14
.qi
= old
[1];
257 newparam
->ljc14
.qj
= old
[2];
258 set_ljparams(comb
,reppow
,old
[3],old
[4],&newparam
->ljc14
.c6
,&newparam
->ljc14
.c12
);
261 newparam
->ljcnb
.qi
= old
[0];
262 newparam
->ljcnb
.qj
= old
[1];
263 set_ljparams(comb
,reppow
,old
[2],old
[3],&newparam
->ljcnb
.c6
,&newparam
->ljcnb
.c12
);
266 set_ljparams(comb
,reppow
,old
[0],old
[1],&newparam
->lj
.c6
,&newparam
->lj
.c12
);
272 newparam
->pdihs
.phiA
= old
[0];
273 newparam
->pdihs
.cpA
= old
[1];
275 /* Change 20100720: Amber occasionally uses negative multiplicities (mathematically OK),
276 * so I have changed the lower limit to -99 /EL
278 newparam
->pdihs
.phiB
= old
[3];
279 newparam
->pdihs
.cpB
= old
[4];
280 /* If both force constants are zero there is no interaction. Return -1 to signal
281 * this entry should NOT be added.
283 if( fabs(newparam
->pdihs
.cpA
) < GMX_REAL_MIN
&& fabs(newparam
->pdihs
.cpB
) < GMX_REAL_MIN
)
288 newparam
->pdihs
.mult
= round_check(old
[2],-99,ftype
,"multiplicity");
292 newparam
->posres
.fcA
[XX
] = old
[0];
293 newparam
->posres
.fcA
[YY
] = old
[1];
294 newparam
->posres
.fcA
[ZZ
] = old
[2];
295 newparam
->posres
.fcB
[XX
] = old
[3];
296 newparam
->posres
.fcB
[YY
] = old
[4];
297 newparam
->posres
.fcB
[ZZ
] = old
[5];
298 newparam
->posres
.pos0A
[XX
] = old
[6];
299 newparam
->posres
.pos0A
[YY
] = old
[7];
300 newparam
->posres
.pos0A
[ZZ
] = old
[8];
301 newparam
->posres
.pos0B
[XX
] = old
[9];
302 newparam
->posres
.pos0B
[YY
] = old
[10];
303 newparam
->posres
.pos0B
[ZZ
] = old
[11];
306 newparam
->disres
.label
= round_check(old
[0],0,ftype
,"label");
307 newparam
->disres
.type
= round_check(old
[1],1,ftype
,"type'");
308 newparam
->disres
.low
= old
[2];
309 newparam
->disres
.up1
= old
[3];
310 newparam
->disres
.up2
= old
[4];
311 newparam
->disres
.kfac
= old
[5];
314 newparam
->orires
.ex
= round_check(old
[0],1,ftype
,"experiment") - 1;
315 newparam
->orires
.label
= round_check(old
[1],1,ftype
,"label");
316 newparam
->orires
.power
= round_check(old
[2],0,ftype
,"power");
317 newparam
->orires
.c
= old
[3];
318 newparam
->orires
.obs
= old
[4];
319 newparam
->orires
.kfac
= old
[5];
322 newparam
->dihres
.phiA
= old
[0];
323 newparam
->dihres
.dphiA
= old
[1];
324 newparam
->dihres
.kfacA
= old
[2];
325 newparam
->dihres
.phiB
= old
[3];
326 newparam
->dihres
.dphiB
= old
[4];
327 newparam
->dihres
.kfacB
= old
[5];
330 for (i
=0; (i
<NR_RBDIHS
); i
++) {
331 newparam
->rbdihs
.rbcA
[i
]=old
[i
];
332 newparam
->rbdihs
.rbcB
[i
]=old
[NR_RBDIHS
+i
];
336 /* Read the dihedral parameters to temporary arrays,
337 * and convert them to the computationally faster
338 * Ryckaert-Bellemans form.
340 /* Use conversion formula for OPLS to Ryckaert-Bellemans: */
341 newparam
->rbdihs
.rbcA
[0]=old
[1]+0.5*(old
[0]+old
[2]);
342 newparam
->rbdihs
.rbcA
[1]=0.5*(3.0*old
[2]-old
[0]);
343 newparam
->rbdihs
.rbcA
[2]=4.0*old
[3]-old
[1];
344 newparam
->rbdihs
.rbcA
[3]=-2.0*old
[2];
345 newparam
->rbdihs
.rbcA
[4]=-4.0*old
[3];
346 newparam
->rbdihs
.rbcA
[5]=0.0;
348 newparam
->rbdihs
.rbcB
[0]=old
[NR_FOURDIHS
+1]+0.5*(old
[NR_FOURDIHS
+0]+old
[NR_FOURDIHS
+2]);
349 newparam
->rbdihs
.rbcB
[1]=0.5*(3.0*old
[NR_FOURDIHS
+2]-old
[NR_FOURDIHS
+0]);
350 newparam
->rbdihs
.rbcB
[2]=4.0*old
[NR_FOURDIHS
+3]-old
[NR_FOURDIHS
+1];
351 newparam
->rbdihs
.rbcB
[3]=-2.0*old
[NR_FOURDIHS
+2];
352 newparam
->rbdihs
.rbcB
[4]=-4.0*old
[NR_FOURDIHS
+3];
353 newparam
->rbdihs
.rbcB
[5]=0.0;
357 newparam
->constr
.dA
= old
[0];
358 newparam
->constr
.dB
= old
[1];
361 newparam
->settle
.doh
=old
[0];
362 newparam
->settle
.dhh
=old
[1];
370 newparam
->vsite
.a
=old
[0];
371 newparam
->vsite
.b
=old
[1];
372 newparam
->vsite
.c
=old
[2];
373 newparam
->vsite
.d
=old
[3];
374 newparam
->vsite
.e
=old
[4];
375 newparam
->vsite
.f
=old
[5];
378 newparam
->vsite
.a
=old
[1] * cos(DEG2RAD
* old
[0]);
379 newparam
->vsite
.b
=old
[1] * sin(DEG2RAD
* old
[0]);
380 newparam
->vsite
.c
=old
[2];
381 newparam
->vsite
.d
=old
[3];
382 newparam
->vsite
.e
=old
[4];
383 newparam
->vsite
.f
=old
[5];
386 newparam
->vsiten
.n
= round_check(old
[0],1,ftype
,"number of atoms");
387 newparam
->vsiten
.a
= old
[1];
390 newparam
->cmap
.cmapA
=old
[0];
391 newparam
->cmap
.cmapB
=old
[1];
396 newparam
->gb
.sar
= old
[0];
397 newparam
->gb
.st
= old
[1];
398 newparam
->gb
.pi
= old
[2];
399 newparam
->gb
.gbr
= old
[3];
400 newparam
->gb
.bmlt
= old
[4];
403 gmx_fatal(FARGS
,"unknown function type %d in %s line %d",
404 ftype
,__FILE__
,__LINE__
);
409 static int enter_params(gmx_ffparams_t
*ffparams
, t_functype ftype
,
410 real forceparams
[MAXFORCEPARAM
],int comb
,real reppow
,
411 int start
,gmx_bool bAppend
)
417 if( (rc
=assign_param(ftype
,&newparam
,forceparams
,comb
,reppow
))<0 )
419 /* -1 means this interaction is all-zero and should not be added */
424 for (type
=start
; (type
<ffparams
->ntypes
); type
++) {
425 if (ffparams
->functype
[type
]==ftype
) {
426 if (F_GB13
== ftype
) {
427 /* Occasionally, the way the 1-3 reference distance is
428 * computed can lead to non-binary-identical results, but I
430 if ((gmx_within_tol(newparam
.gb
.sar
, ffparams
->iparams
[type
].gb
.sar
, 1e-6)) &&
431 (gmx_within_tol(newparam
.gb
.st
, ffparams
->iparams
[type
].gb
.st
, 1e-6)) &&
432 (gmx_within_tol(newparam
.gb
.pi
, ffparams
->iparams
[type
].gb
.pi
, 1e-6)) &&
433 (gmx_within_tol(newparam
.gb
.gbr
, ffparams
->iparams
[type
].gb
.gbr
, 1e-6)) &&
434 (gmx_within_tol(newparam
.gb
.bmlt
, ffparams
->iparams
[type
].gb
.bmlt
, 1e-6))) {
439 if (memcmp(&newparam
,&ffparams
->iparams
[type
],(size_t)sizeof(newparam
)) == 0)
446 type
= ffparams
->ntypes
;
449 fprintf(debug
,"copying newparam to ffparams->iparams[%d] (ntypes=%d)\n",
450 type
,ffparams
->ntypes
);
451 memcpy(&ffparams
->iparams
[type
],&newparam
,(size_t)sizeof(newparam
));
454 ffparams
->functype
[type
]=ftype
;
459 static void append_interaction(t_ilist
*ilist
,
460 int type
,int nral
,atom_id a
[MAXATOMLIST
])
467 ilist
->iatoms
[where1
++]=type
;
468 for (i
=0; (i
<nral
); i
++)
469 ilist
->iatoms
[where1
++]=a
[i
];
472 static void enter_function(t_params
*p
,t_functype ftype
,int comb
,real reppow
,
473 gmx_ffparams_t
*ffparams
,t_ilist
*il
,
475 gmx_bool bNB
,gmx_bool bAppend
)
477 int k
,type
,nr
,nral
,delta
,start
;
479 start
= ffparams
->ntypes
;
482 for (k
=0; k
<nr
; k
++) {
483 if (*maxtypes
<= ffparams
->ntypes
) {
485 srenew(ffparams
->functype
,*maxtypes
);
486 srenew(ffparams
->iparams
, *maxtypes
);
488 fprintf(debug
,"%s, line %d: srenewed idef->functype and idef->iparams to %d\n",
489 __FILE__
,__LINE__
,*maxtypes
);
491 type
= enter_params(ffparams
,ftype
,p
->param
[k
].c
,comb
,reppow
,start
,bAppend
);
492 /* Type==-1 is used as a signal that this interaction is all-zero and should not be added. */
493 if (!bNB
&& type
>=0) {
496 srenew(il
->iatoms
,il
->nr
+delta
);
497 append_interaction(il
,type
,nral
,p
->param
[k
].a
);
502 void convert_params(int atnr
,t_params nbtypes
[],
503 t_molinfo
*mi
,int comb
,double reppow
,real fudgeQQ
,
514 ffp
= &mtop
->ffparams
;
517 ffp
->functype
= NULL
;
519 ffp
->reppow
= reppow
;
521 enter_function(&(nbtypes
[F_LJ
]), (t_functype
)F_LJ
, comb
,reppow
,ffp
,NULL
,
522 &maxtypes
,TRUE
,TRUE
);
523 enter_function(&(nbtypes
[F_BHAM
]),(t_functype
)F_BHAM
, comb
,reppow
,ffp
,NULL
,
524 &maxtypes
,TRUE
,TRUE
);
526 for(mt
=0; mt
<mtop
->nmoltype
; mt
++) {
527 molt
= &mtop
->moltype
[mt
];
528 for(i
=0; (i
<F_NRE
); i
++) {
529 molt
->ilist
[i
].nr
= 0;
530 molt
->ilist
[i
].iatoms
= NULL
;
532 plist
= mi
[mt
].plist
;
534 flags
= interaction_function
[i
].flags
;
535 if ((i
!= F_LJ
) && (i
!= F_BHAM
) && ((flags
& IF_BOND
) ||
536 (flags
& IF_VSITE
) ||
537 (flags
& IF_CONSTRAINT
))) {
538 enter_function(&(plist
[i
]),(t_functype
)i
,comb
,reppow
,
540 &maxtypes
,FALSE
,(i
== F_POSRES
));
545 fprintf(debug
,"%s, line %d: There are %d functypes in idef\n",
546 __FILE__
,__LINE__
,ffp
->ntypes
);
549 ffp
->fudgeQQ
= fudgeQQ
;