3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
35 /* This file is completely threadsafe - keep it that way! */
46 #include "gmx_fatal.h"
51 #include "gpp_atomtype.h"
53 static int round_check(real r
,int limit
,int ftype
,const char *name
)
62 if (r
-i
> 0.01 || r
-i
< -0.01)
63 gmx_fatal(FARGS
,"A non-integer value (%f) was supplied for '%s' in %s",
64 r
,name
,interaction_function
[ftype
].longname
);
67 gmx_fatal(FARGS
,"Value of '%s' in %s is %d, which is smaller than the minimum of %d",
68 name
,interaction_function
[ftype
].longname
,i
,limit
);
73 static void set_ljparams(int comb
,double reppow
,real v
,real w
,
76 if (comb
== eCOMB_ARITHMETIC
|| comb
== eCOMB_GEOM_SIG_EPS
) {
79 *c12
= 4*w
*pow(v
,reppow
);
81 /* Interpret negative sigma as c6=0 and c12 with -sigma */
83 *c12
= 4*w
*pow(-v
,reppow
);
91 static void assign_param(t_functype ftype
,t_iparams
*newparam
,
92 real old
[MAXFORCEPARAM
],int comb
,double reppow
)
98 for(j
=0; (j
<MAXFORCEPARAM
); j
++)
100 newparam
->generic
.buf
[j
]=0.0;
104 /* Post processing of input data: store cosine iso angle itself */
105 newparam
->harmonic
.rA
=cos(old
[0]*DEG2RAD
);
106 newparam
->harmonic
.krA
=old
[1];
107 newparam
->harmonic
.rB
=cos(old
[2]*DEG2RAD
);
108 newparam
->harmonic
.krB
=old
[3];
111 /* Post processing of input data: store square of length itself */
112 newparam
->harmonic
.rA
=sqr(old
[0]);
113 newparam
->harmonic
.krA
=old
[1];
114 newparam
->harmonic
.rB
=sqr(old
[2]);
115 newparam
->harmonic
.krB
=old
[3];
118 newparam
->fene
.bm
=old
[0];
119 newparam
->fene
.kb
=old
[1];
122 newparam
->restraint
.lowA
= old
[0];
123 newparam
->restraint
.up1A
= old
[1];
124 newparam
->restraint
.up2A
= old
[2];
125 newparam
->restraint
.kA
= old
[3];
126 newparam
->restraint
.lowB
= old
[4];
127 newparam
->restraint
.up1B
= old
[5];
128 newparam
->restraint
.up2B
= old
[6];
129 newparam
->restraint
.kB
= old
[7];
135 newparam
->tab
.table
= round_check(old
[0],0,ftype
,"table index");
136 newparam
->tab
.kA
= old
[1];
137 newparam
->tab
.kB
= old
[3];
139 case F_CROSS_BOND_BONDS
:
140 newparam
->cross_bb
.r1e
=old
[0];
141 newparam
->cross_bb
.r2e
=old
[1];
142 newparam
->cross_bb
.krr
=old
[2];
144 case F_CROSS_BOND_ANGLES
:
145 newparam
->cross_ba
.r1e
=old
[0];
146 newparam
->cross_ba
.r2e
=old
[1];
147 newparam
->cross_ba
.r3e
=old
[2];
148 newparam
->cross_ba
.krt
=old
[3];
151 newparam
->u_b
.theta
=old
[0];
152 newparam
->u_b
.ktheta
=old
[1];
153 newparam
->u_b
.r13
=old
[2];
154 newparam
->u_b
.kUB
=old
[3];
156 case F_QUARTIC_ANGLES
:
157 newparam
->qangle
.theta
=old
[0];
159 newparam
->qangle
.c
[i
]=old
[i
+1];
165 newparam
->harmonic
.rA
=old
[0];
166 newparam
->harmonic
.krA
=old
[1];
167 newparam
->harmonic
.rB
=old
[2];
168 newparam
->harmonic
.krB
=old
[3];
171 newparam
->morse
.b0
=old
[0];
172 newparam
->morse
.cb
=old
[1];
173 newparam
->morse
.beta
=old
[2];
176 newparam
->cubic
.b0
=old
[0];
177 newparam
->cubic
.kb
=old
[1];
178 newparam
->cubic
.kcub
=old
[2];
183 newparam
->polarize
.alpha
= old
[0];
186 newparam
->wpol
.al_x
=old
[0];
187 newparam
->wpol
.al_y
=old
[1];
188 newparam
->wpol
.al_z
=old
[2];
189 newparam
->wpol
.rOH
=old
[3];
190 newparam
->wpol
.rHH
=old
[4];
191 newparam
->wpol
.rOD
=old
[5];
194 newparam
->thole
.a
= old
[0];
195 newparam
->thole
.alpha1
= old
[1];
196 newparam
->thole
.alpha2
= old
[2];
197 if ((old
[1] > 0) && (old
[2] > 0))
198 newparam
->thole
.rfac
= old
[0]*pow(old
[1]*old
[2],-1.0/6.0);
200 newparam
->thole
.rfac
= 1;
203 newparam
->bham
.a
= old
[0];
204 newparam
->bham
.b
= old
[1];
205 newparam
->bham
.c
= old
[2];
208 set_ljparams(comb
,reppow
,old
[0],old
[1],&newparam
->lj14
.c6A
,&newparam
->lj14
.c12A
);
209 set_ljparams(comb
,reppow
,old
[2],old
[3],&newparam
->lj14
.c6B
,&newparam
->lj14
.c12B
);
212 newparam
->ljc14
.fqq
= old
[0];
213 newparam
->ljc14
.qi
= old
[1];
214 newparam
->ljc14
.qj
= old
[2];
215 set_ljparams(comb
,reppow
,old
[3],old
[4],&newparam
->ljc14
.c6
,&newparam
->ljc14
.c12
);
218 newparam
->ljcnb
.qi
= old
[0];
219 newparam
->ljcnb
.qj
= old
[1];
220 set_ljparams(comb
,reppow
,old
[2],old
[3],&newparam
->ljcnb
.c6
,&newparam
->ljcnb
.c12
);
223 set_ljparams(comb
,reppow
,old
[0],old
[1],&newparam
->lj
.c6
,&newparam
->lj
.c12
);
228 newparam
->pdihs
.phiA
= old
[0];
229 newparam
->pdihs
.cpA
= old
[1];
231 /* Dont do any checks if all parameters are zero (such interactions will be removed) */
232 tmp
=fabs(old
[0])+fabs(old
[1])+fabs(old
[2])+fabs(old
[3])+fabs(old
[4]);
233 newparam
->pdihs
.mult
= (tmp
< GMX_REAL_MIN
) ? 0 : round_check(old
[2],1,ftype
,"multiplicity");
235 newparam
->pdihs
.phiB
= old
[3];
236 newparam
->pdihs
.cpB
= old
[4];
239 newparam
->posres
.fcA
[XX
] = old
[0];
240 newparam
->posres
.fcA
[YY
] = old
[1];
241 newparam
->posres
.fcA
[ZZ
] = old
[2];
242 newparam
->posres
.fcB
[XX
] = old
[3];
243 newparam
->posres
.fcB
[YY
] = old
[4];
244 newparam
->posres
.fcB
[ZZ
] = old
[5];
245 newparam
->posres
.pos0A
[XX
] = old
[6];
246 newparam
->posres
.pos0A
[YY
] = old
[7];
247 newparam
->posres
.pos0A
[ZZ
] = old
[8];
248 newparam
->posres
.pos0B
[XX
] = old
[9];
249 newparam
->posres
.pos0B
[YY
] = old
[10];
250 newparam
->posres
.pos0B
[ZZ
] = old
[11];
253 newparam
->disres
.label
= round_check(old
[0],0,ftype
,"label");
254 newparam
->disres
.type
= round_check(old
[1],1,ftype
,"type'");
255 newparam
->disres
.low
= old
[2];
256 newparam
->disres
.up1
= old
[3];
257 newparam
->disres
.up2
= old
[4];
258 newparam
->disres
.kfac
= old
[5];
261 newparam
->orires
.ex
= round_check(old
[0],1,ftype
,"experiment") - 1;
262 newparam
->orires
.label
= round_check(old
[1],1,ftype
,"label");
263 newparam
->orires
.power
= round_check(old
[2],0,ftype
,"power");
264 newparam
->orires
.c
= old
[3];
265 newparam
->orires
.obs
= old
[4];
266 newparam
->orires
.kfac
= old
[5];
269 newparam
->dihres
.label
= round_check(old
[0],0,ftype
,"label");
270 newparam
->dihres
.phi
= old
[1];
271 newparam
->dihres
.dphi
= old
[2];
272 newparam
->dihres
.kfac
= old
[3];
273 newparam
->dihres
.power
= round_check(old
[4],0,ftype
,"power");
276 for (i
=0; (i
<NR_RBDIHS
); i
++) {
277 newparam
->rbdihs
.rbcA
[i
]=old
[i
];
278 newparam
->rbdihs
.rbcB
[i
]=old
[NR_RBDIHS
+i
];
282 /* Read the dihedral parameters to temporary arrays,
283 * and convert them to the computationally faster
284 * Ryckaert-Bellemans form.
286 /* Use conversion formula for OPLS to Ryckaert-Bellemans: */
287 newparam
->rbdihs
.rbcA
[0]=old
[1]+0.5*(old
[0]+old
[2]);
288 newparam
->rbdihs
.rbcA
[1]=0.5*(3.0*old
[2]-old
[0]);
289 newparam
->rbdihs
.rbcA
[2]=4.0*old
[3]-old
[1];
290 newparam
->rbdihs
.rbcA
[3]=-2.0*old
[2];
291 newparam
->rbdihs
.rbcA
[4]=-4.0*old
[3];
292 newparam
->rbdihs
.rbcA
[5]=0.0;
294 newparam
->rbdihs
.rbcB
[0]=old
[NR_FOURDIHS
+1]+0.5*(old
[NR_FOURDIHS
+0]+old
[NR_FOURDIHS
+2]);
295 newparam
->rbdihs
.rbcB
[1]=0.5*(3.0*old
[NR_FOURDIHS
+2]-old
[NR_FOURDIHS
+0]);
296 newparam
->rbdihs
.rbcB
[2]=4.0*old
[NR_FOURDIHS
+3]-old
[NR_FOURDIHS
+1];
297 newparam
->rbdihs
.rbcB
[3]=-2.0*old
[NR_FOURDIHS
+2];
298 newparam
->rbdihs
.rbcB
[4]=-4.0*old
[NR_FOURDIHS
+3];
299 newparam
->rbdihs
.rbcB
[5]=0.0;
303 newparam
->constr
.dA
= old
[0];
304 newparam
->constr
.dB
= old
[1];
307 newparam
->settle
.doh
=old
[0];
308 newparam
->settle
.dhh
=old
[1];
316 newparam
->vsite
.a
=old
[0];
317 newparam
->vsite
.b
=old
[1];
318 newparam
->vsite
.c
=old
[2];
319 newparam
->vsite
.d
=old
[3];
320 newparam
->vsite
.e
=old
[4];
321 newparam
->vsite
.f
=old
[5];
324 newparam
->vsite
.a
=old
[1] * cos(DEG2RAD
* old
[0]);
325 newparam
->vsite
.b
=old
[1] * sin(DEG2RAD
* old
[0]);
326 newparam
->vsite
.c
=old
[2];
327 newparam
->vsite
.d
=old
[3];
328 newparam
->vsite
.e
=old
[4];
329 newparam
->vsite
.f
=old
[5];
332 newparam
->vsiten
.n
= round_check(old
[0],1,ftype
,"number of atoms");
333 newparam
->vsiten
.a
= old
[1];
336 newparam
->cmap
.cmapA
=old
[0];
337 newparam
->cmap
.cmapB
=old
[1];
342 newparam
->gb
.sar
= old
[0];
343 newparam
->gb
.st
= old
[1];
344 newparam
->gb
.pi
= old
[2];
345 newparam
->gb
.gbr
= old
[3];
346 newparam
->gb
.bmlt
= old
[4];
349 gmx_fatal(FARGS
,"unknown function type %d in %s line %d",
350 ftype
,__FILE__
,__LINE__
);
354 static int enter_params(gmx_ffparams_t
*ffparams
, t_functype ftype
,
355 real forceparams
[MAXFORCEPARAM
],int comb
,real reppow
,
356 int start
,bool bAppend
)
361 assign_param(ftype
,&newparam
,forceparams
,comb
,reppow
);
363 for (type
=start
; (type
<ffparams
->ntypes
); type
++) {
364 if (ffparams
->functype
[type
]==ftype
) {
365 if (memcmp(&newparam
,&ffparams
->iparams
[type
],(size_t)sizeof(newparam
)) == 0)
371 type
= ffparams
->ntypes
;
374 fprintf(debug
,"copying newparam to ffparams->iparams[%d] (ntypes=%d)\n",
375 type
,ffparams
->ntypes
);
376 memcpy(&ffparams
->iparams
[type
],&newparam
,(size_t)sizeof(newparam
));
379 ffparams
->functype
[type
]=ftype
;
384 static void append_interaction(t_ilist
*ilist
,
385 int type
,int nral
,atom_id a
[MAXATOMLIST
])
392 ilist
->iatoms
[where1
++]=type
;
393 for (i
=0; (i
<nral
); i
++)
394 ilist
->iatoms
[where1
++]=a
[i
];
397 static void enter_function(t_params
*p
,t_functype ftype
,int comb
,real reppow
,
398 gmx_ffparams_t
*ffparams
,t_ilist
*il
,
400 bool bNB
,bool bAppend
)
402 int k
,type
,nr
,nral
,delta
,start
;
404 start
= ffparams
->ntypes
;
407 for (k
=0; k
<nr
; k
++) {
408 if (*maxtypes
<= ffparams
->ntypes
) {
410 srenew(ffparams
->functype
,*maxtypes
);
411 srenew(ffparams
->iparams
, *maxtypes
);
413 fprintf(debug
,"%s, line %d: srenewed idef->functype and idef->iparams to %d\n",
414 __FILE__
,__LINE__
,*maxtypes
);
416 type
= enter_params(ffparams
,ftype
,p
->param
[k
].c
,comb
,reppow
,start
,bAppend
);
420 srenew(il
->iatoms
,il
->nr
+delta
);
421 append_interaction(il
,type
,nral
,p
->param
[k
].a
);
426 static void new_interaction_list(t_ilist
*ilist
)
434 void convert_params(int atnr
,t_params nbtypes
[],
435 t_molinfo
*mi
,int comb
,double reppow
,real fudgeQQ
,
446 ffp
= &mtop
->ffparams
;
449 ffp
->functype
= NULL
;
451 ffp
->reppow
= reppow
;
453 enter_function(&(nbtypes
[F_LJ
]), (t_functype
)F_LJ
, comb
,reppow
,ffp
,NULL
,
454 &maxtypes
,TRUE
,TRUE
);
455 enter_function(&(nbtypes
[F_BHAM
]),(t_functype
)F_BHAM
, comb
,reppow
,ffp
,NULL
,
456 &maxtypes
,TRUE
,TRUE
);
458 for(mt
=0; mt
<mtop
->nmoltype
; mt
++) {
459 molt
= &mtop
->moltype
[mt
];
460 for(i
=0; (i
<F_NRE
); i
++) {
461 molt
->ilist
[i
].nr
= 0;
462 molt
->ilist
[i
].iatoms
= NULL
;
464 plist
= mi
[mt
].plist
;
466 flags
= interaction_function
[i
].flags
;
467 if ((i
!= F_LJ
) && (i
!= F_BHAM
) && ((flags
& IF_BOND
) ||
468 (flags
& IF_VSITE
) ||
469 (flags
& IF_CONSTRAINT
))) {
470 enter_function(&(plist
[i
]),(t_functype
)i
,comb
,reppow
,
472 &maxtypes
,FALSE
,(i
== F_POSRES
));
477 fprintf(debug
,"%s, line %d: There are %d functypes in idef\n",
478 __FILE__
,__LINE__
,ffp
->ntypes
);
481 ffp
->fudgeQQ
= fudgeQQ
;