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,2013, 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.
46 #include "vsite_parm.h"
55 #include "gmx_fatal.h"
72 vsitebondparam_t *vsbp;
80 static int vsite_bond_nrcheck(int ftype)
84 if ((interaction_function[ftype].flags & (IF_BTYPE | IF_CONSTRAINT | IF_ATYPE)) || (ftype == F_IDIHS))
85 nrcheck = NRAL(ftype);
92 static void enter_bonded(int nratoms, int *nrbonded, t_mybonded **bondeds,
97 srenew(*bondeds, *nrbonded+1);
99 /* copy atom numbers */
100 for(j=0; j<nratoms; j++)
101 (*bondeds)[*nrbonded].a[j] = param->a[j];
103 (*bondeds)[*nrbonded].c = param->C0;
108 static void get_bondeds(int nrat, t_iatom atoms[],
109 at2vsitebond_t *at2vb, t_params plist[],
110 int *nrbond, t_mybonded **bonds,
111 int *nrang, t_mybonded **angles,
112 int *nridih, t_mybonded **idihs )
114 int k,i,ftype,nrcheck;
117 for(k=0; k<nrat; k++) {
118 for(i=0; i<at2vb[atoms[k]].nr; i++) {
119 ftype = at2vb[atoms[k]].vsbp[i].ftype;
120 param = at2vb[atoms[k]].vsbp[i].param;
121 nrcheck = vsite_bond_nrcheck(ftype);
122 /* abuse nrcheck to see if we're adding bond, angle or idih */
124 case 2: enter_bonded(nrcheck,nrbond,bonds, param); break;
125 case 3: enter_bonded(nrcheck,nrang, angles,param); break;
126 case 4: enter_bonded(nrcheck,nridih,idihs, param); break;
132 static at2vsitebond_t *make_at2vsitebond(int natoms,t_params plist[])
135 int ftype,i,j,nrcheck,nr;
137 at2vsitebond_t *at2vb;
142 for(ftype=0; (ftype<F_NRE); ftype++) {
143 if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN) {
144 for(i=0; (i<plist[ftype].nr); i++) {
145 for(j=0; j<NRAL(ftype); j++)
146 bVSI[plist[ftype].param[i].a[j]] = TRUE;
151 for(ftype=0; (ftype<F_NRE); ftype++) {
152 nrcheck = vsite_bond_nrcheck(ftype);
154 for(i=0; (i<plist[ftype].nr); i++) {
155 aa = plist[ftype].param[i].a;
156 for(j=0; j<nrcheck; j++) {
158 nr = at2vb[aa[j]].nr;
160 srenew(at2vb[aa[j]].vsbp,nr+10);
161 at2vb[aa[j]].vsbp[nr].ftype = ftype;
162 at2vb[aa[j]].vsbp[nr].param = &plist[ftype].param[i];
174 static void done_at2vsitebond(int natoms,at2vsitebond_t *at2vb)
178 for(i=0; i<natoms; i++)
180 sfree(at2vb[i].vsbp);
184 static at2vsitecon_t *make_at2vsitecon(int natoms,t_params plist[])
187 int ftype,i,j,ai,aj,nr;
188 at2vsitecon_t *at2vc;
193 for(ftype=0; (ftype<F_NRE); ftype++) {
194 if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN) {
195 for(i=0; (i<plist[ftype].nr); i++) {
196 for(j=0; j<NRAL(ftype); j++)
197 bVSI[plist[ftype].param[i].a[j]] = TRUE;
202 for(ftype=0; (ftype<F_NRE); ftype++) {
203 if (interaction_function[ftype].flags & IF_CONSTRAINT) {
204 for(i=0; (i<plist[ftype].nr); i++) {
205 ai = plist[ftype].param[i].AI;
206 aj = plist[ftype].param[i].AJ;
207 if (bVSI[ai] && bVSI[aj]) {
208 /* Store forward direction */
211 srenew(at2vc[ai].aj,nr+10);
212 at2vc[ai].aj[nr] = aj;
214 /* Store backward direction */
217 srenew(at2vc[aj].aj,nr+10);
218 at2vc[aj].aj[nr] = ai;
229 static void done_at2vsitecon(int natoms,at2vsitecon_t *at2vc)
233 for(i=0; i<natoms; i++)
240 static void print_bad(FILE *fp,
241 int nrbond, t_mybonded *bonds,
242 int nrang, t_mybonded *angles,
243 int nridih, t_mybonded *idihs )
248 fprintf(fp,"bonds:");
249 for(i=0; i<nrbond; i++)
250 fprintf(fp," %u-%u (%g)",
251 bonds[i].AI+1, bonds[i].AJ+1, bonds[i].c);
255 fprintf(fp,"angles:");
256 for(i=0; i<nrang; i++)
257 fprintf(fp," %u-%u-%u (%g)",
258 angles[i].AI+1, angles[i].AJ+1,
259 angles[i].AK+1, angles[i].c);
263 fprintf(fp,"idihs:");
264 for(i=0; i<nridih; i++)
265 fprintf(fp," %u-%u-%u-%u (%g)",
266 idihs[i].AI+1, idihs[i].AJ+1,
267 idihs[i].AK+1, idihs[i].AL+1, idihs[i].c);
272 static void print_param(FILE *fp, int ftype, int i, t_param *param)
275 static int prev_ftype= NOTSET;
276 static int prev_i = NOTSET;
279 if ( (ftype!=prev_ftype) || (i!=prev_i) ) {
284 fprintf(fp,"(%d) plist[%s].param[%d]",
285 pass,interaction_function[ftype].name,i);
286 for(j=0; j<NRFP(ftype); j++)
287 fprintf(fp,".c[%d]=%g ",j,param->c[j]);
292 static real get_bond_length(int nrbond, t_mybonded bonds[],
293 t_iatom ai, t_iatom aj)
299 for (i=0; i < nrbond && (bondlen==NOTSET); i++) {
300 /* check both ways */
301 if ( ( (ai == bonds[i].AI) && (aj == bonds[i].AJ) ) ||
302 ( (ai == bonds[i].AJ) && (aj == bonds[i].AI) ) )
303 bondlen = bonds[i].c; /* note: bonds[i].c might be NOTSET */
308 static real get_angle(int nrang, t_mybonded angles[],
309 t_iatom ai, t_iatom aj, t_iatom ak)
315 for (i=0; i < nrang && (angle==NOTSET); i++) {
316 /* check both ways */
317 if ( ( (ai==angles[i].AI) && (aj==angles[i].AJ) && (ak==angles[i].AK) ) ||
318 ( (ai==angles[i].AK) && (aj==angles[i].AJ) && (ak==angles[i].AI) ) )
319 angle = DEG2RAD*angles[i].c;
324 static char *get_atomtype_name_AB(t_atom *atom,gpp_atomtype_t atype)
328 name = get_atomtype_name(atom->type,atype);
330 /* When using the decoupling option, atom types are changed
331 * to decoupled for the non-bonded interactions, but the virtual
332 * sites constructions should be based on the "normal" interactions.
333 * So we return the state B atom type name if the state A atom
334 * type is the decoupled one. We should actually check for the atom
335 * type number, but that's not passed here. So we check for
336 * the decoupled atom type name. This should not cause trouble
337 * as this code is only used for topologies with v-sites without
338 * parameters generated by pdb2gmx.
340 if (strcmp(name,"decoupled") == 0)
342 name = get_atomtype_name(atom->typeB,atype);
348 static gmx_bool calc_vsite3_param(gpp_atomtype_t atype,
349 t_param *param, t_atoms *at,
350 int nrbond, t_mybonded *bonds,
351 int nrang, t_mybonded *angles )
353 /* i = virtual site | ,k
354 * j = 1st bonded heavy atom | i-j
355 * k,l = 2nd bonded atoms | `l
358 gmx_bool bXH3,bError;
359 real bjk,bjl,a=-1,b=-1;
360 /* check if this is part of a NH3 , NH2-umbrella or CH3 group,
361 * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
365 fprintf(debug,"atom %u type %s ",
367 get_atomtype_name_AB(&at->atom[param->a[i]],atype));
371 ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AK],atype),"MNH",3)==0) &&
372 (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AL],atype),"MNH",3)==0) ) ||
373 ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AK],atype),"MCH3",4)==0) &&
374 (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AL],atype),"MCH3",4)==0) );
376 bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
377 bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
378 bError = (bjk==NOTSET) || (bjl==NOTSET);
380 /* now we get some XH2/XH3 group specific construction */
381 /* note: we call the heavy atom 'C' and the X atom 'N' */
382 real bMM,bCM,bCN,bNH,aCNH,dH,rH,dM,rM;
385 /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
386 bError = bError || (bjk!=bjl);
388 /* the X atom (C or N) in the XH2/XH3 group is the first after the masses: */
389 aN = max(param->AK,param->AL)+1;
391 /* get common bonds */
392 bMM = get_bond_length(nrbond, bonds, param->AK, param->AL);
394 bCN = get_bond_length(nrbond, bonds, param->AJ, aN);
395 bError = bError || (bMM==NOTSET) || (bCN==NOTSET);
397 /* calculate common things */
399 dM = sqrt( sqr(bCM) - sqr(rM) );
401 /* are we dealing with the X atom? */
402 if ( param->AI == aN ) {
403 /* this is trivial */
404 a = b = 0.5 * bCN/dM;
407 /* get other bondlengths and angles: */
408 bNH = get_bond_length(nrbond, bonds, aN, param->AI);
409 aCNH= get_angle (nrang, angles, param->AJ, aN, param->AI);
410 bError = bError || (bNH==NOTSET) || (aCNH==NOTSET);
413 dH = bCN - bNH * cos(aCNH);
414 rH = bNH * sin(aCNH);
416 a = 0.5 * ( dH/dM + rH/rM );
417 b = 0.5 * ( dH/dM - rH/rM );
420 gmx_fatal(FARGS,"calc_vsite3_param not implemented for the general case "
421 "(atom %d)",param->AI+1);
427 fprintf(debug,"params for vsite3 %u: %g %g\n",
428 param->AI+1,param->C0,param->C1);
433 static gmx_bool calc_vsite3fd_param(t_param *param,
434 int nrbond, t_mybonded *bonds,
435 int nrang, t_mybonded *angles)
437 /* i = virtual site | ,k
438 * j = 1st bonded heavy atom | i-j
439 * k,l = 2nd bonded atoms | `l
443 real bij,bjk,bjl,aijk,aijl,rk,rl;
445 bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
446 bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
447 bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
448 aijk= get_angle (nrang, angles, param->AI, param->AJ, param->AK);
449 aijl= get_angle (nrang, angles, param->AI, param->AJ, param->AL);
450 bError = (bij==NOTSET) || (bjk==NOTSET) || (bjl==NOTSET) ||
451 (aijk==NOTSET) || (aijl==NOTSET);
453 rk = bjk * sin(aijk);
454 rl = bjl * sin(aijl);
455 param->C0 = rk / (rk + rl);
456 param->C1 = -bij; /* 'bond'-length for fixed distance vsite */
459 fprintf(debug,"params for vsite3fd %u: %g %g\n",
460 param->AI+1,param->C0,param->C1);
464 static gmx_bool calc_vsite3fad_param(t_param *param,
465 int nrbond, t_mybonded *bonds,
466 int nrang, t_mybonded *angles)
468 /* i = virtual site |
469 * j = 1st bonded heavy atom | i-j
470 * k = 2nd bonded heavy atom | `k-l
471 * l = 3d bonded heavy atom |
474 gmx_bool bSwapParity,bError;
477 bSwapParity = ( param->C1 == -1 );
479 bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
480 aijk = get_angle (nrang, angles, param->AI, param->AJ, param->AK);
481 bError = (bij==NOTSET) || (aijk==NOTSET);
483 param->C1 = bij; /* 'bond'-length for fixed distance vsite */
484 param->C0 = RAD2DEG*aijk; /* 'bond'-angle for fixed angle vsite */
487 param->C0 = 360 - param->C0;
490 fprintf(debug,"params for vsite3fad %u: %g %g\n",
491 param->AI+1,param->C0,param->C1);
495 static gmx_bool calc_vsite3out_param(gpp_atomtype_t atype,
496 t_param *param, t_atoms *at,
497 int nrbond, t_mybonded *bonds,
498 int nrang, t_mybonded *angles)
500 /* i = virtual site | ,k
501 * j = 1st bonded heavy atom | i-j
502 * k,l = 2nd bonded atoms | `l
503 * NOTE: i is out of the j-k-l plane!
506 gmx_bool bXH3,bError,bSwapParity;
507 real bij,bjk,bjl,aijk,aijl,akjl,pijk,pijl,a,b,c;
509 /* check if this is part of a NH2-umbrella, NH3 or CH3 group,
510 * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
514 fprintf(debug,"atom %u type %s ",
515 param->a[i]+1,get_atomtype_name_AB(&at->atom[param->a[i]],atype));
519 ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AK],atype),"MNH",3)==0) &&
520 (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AL],atype),"MNH",3)==0) ) ||
521 ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AK],atype),"MCH3",4)==0) &&
522 (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AL],atype),"MCH3",4)==0) );
524 /* check if construction parity must be swapped */
525 bSwapParity = ( param->C1 == -1 );
527 bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
528 bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
529 bError = (bjk==NOTSET) || (bjl==NOTSET);
531 /* now we get some XH3 group specific construction */
532 /* note: we call the heavy atom 'C' and the X atom 'N' */
533 real bMM,bCM,bCN,bNH,aCNH,dH,rH,rHx,rHy,dM,rM;
536 /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
537 bError = bError || (bjk!=bjl);
539 /* the X atom (C or N) in the XH3 group is the first after the masses: */
540 aN = max(param->AK,param->AL)+1;
542 /* get all bondlengths and angles: */
543 bMM = get_bond_length(nrbond, bonds, param->AK, param->AL);
545 bCN = get_bond_length(nrbond, bonds, param->AJ, aN);
546 bNH = get_bond_length(nrbond, bonds, aN, param->AI);
547 aCNH= get_angle (nrang, angles, param->AJ, aN, param->AI);
549 (bMM==NOTSET) || (bCN==NOTSET) || (bNH==NOTSET) || (aCNH==NOTSET);
552 dH = bCN - bNH * cos(aCNH);
553 rH = bNH * sin(aCNH);
554 /* we assume the H's are symmetrically distributed */
555 rHx = rH*cos(DEG2RAD*30);
556 rHy = rH*sin(DEG2RAD*30);
558 dM = sqrt( sqr(bCM) - sqr(rM) );
559 a = 0.5*( (dH/dM) - (rHy/rM) );
560 b = 0.5*( (dH/dM) + (rHy/rM) );
564 /* this is the general construction */
566 bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
567 aijk= get_angle (nrang, angles, param->AI, param->AJ, param->AK);
568 aijl= get_angle (nrang, angles, param->AI, param->AJ, param->AL);
569 akjl= get_angle (nrang, angles, param->AK, param->AJ, param->AL);
571 (bij==NOTSET) || (aijk==NOTSET) || (aijl==NOTSET) || (akjl==NOTSET);
573 pijk = cos(aijk)*bij;
574 pijl = cos(aijl)*bij;
575 a = ( pijk + (pijk*cos(akjl)-pijl) * cos(akjl) / sqr(sin(akjl)) ) / bjk;
576 b = ( pijl + (pijl*cos(akjl)-pijk) * cos(akjl) / sqr(sin(akjl)) ) / bjl;
577 c = - sqrt( sqr(bij) -
578 ( sqr(pijk) - 2*pijk*pijl*cos(akjl) + sqr(pijl) )
580 / ( bjk*bjl*sin(akjl) );
590 fprintf(debug,"params for vsite3out %u: %g %g %g\n",
591 param->AI+1,param->C0,param->C1,param->C2);
595 static gmx_bool calc_vsite4fd_param(t_param *param,
596 int nrbond, t_mybonded *bonds,
597 int nrang, t_mybonded *angles)
599 /* i = virtual site | ,k
600 * j = 1st bonded heavy atom | i-j-m
601 * k,l,m = 2nd bonded atoms | `l
605 real bij,bjk,bjl,bjm,aijk,aijl,aijm,akjm,akjl;
606 real pk,pl,pm,cosakl,cosakm,sinakl,sinakm,cl,cm;
608 bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
609 bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
610 bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
611 bjm = get_bond_length(nrbond, bonds, param->AJ, param->AM);
612 aijk= get_angle (nrang, angles, param->AI, param->AJ, param->AK);
613 aijl= get_angle (nrang, angles, param->AI, param->AJ, param->AL);
614 aijm= get_angle (nrang, angles, param->AI, param->AJ, param->AM);
615 akjm= get_angle (nrang, angles, param->AK, param->AJ, param->AM);
616 akjl= get_angle (nrang, angles, param->AK, param->AJ, param->AL);
617 bError = (bij==NOTSET) || (bjk==NOTSET) || (bjl==NOTSET) || (bjm==NOTSET) ||
618 (aijk==NOTSET) || (aijl==NOTSET) || (aijm==NOTSET) || (akjm==NOTSET) ||
625 cosakl = (cos(akjl) - cos(aijk)*cos(aijl)) / (sin(aijk)*sin(aijl));
626 cosakm = (cos(akjm) - cos(aijk)*cos(aijm)) / (sin(aijk)*sin(aijm));
627 if ( cosakl < -1 || cosakl > 1 || cosakm < -1 || cosakm > 1 ) {
628 fprintf(stderr,"virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
629 param->AI+1,RAD2DEG*aijk,RAD2DEG*aijl,RAD2DEG*aijm);
630 gmx_fatal(FARGS,"invalid construction in calc_vsite4fd for atom %d: "
631 "cosakl=%f, cosakm=%f\n",param->AI+1,cosakl,cosakm);
633 sinakl = sqrt(1-sqr(cosakl));
634 sinakm = sqrt(1-sqr(cosakm));
636 /* note: there is a '+' because of the way the sines are calculated */
637 cl = -pk / ( pl*cosakl - pk + pl*sinakl*(pm*cosakm-pk)/(pm*sinakm) );
638 cm = -pk / ( pm*cosakm - pk + pm*sinakm*(pl*cosakl-pk)/(pl*sinakl) );
644 fprintf(debug,"params for vsite4fd %u: %g %g %g\n",
645 param->AI+1,param->C0,param->C1,param->C2);
653 calc_vsite4fdn_param(t_param *param,
654 int nrbond, t_mybonded *bonds,
655 int nrang, t_mybonded *angles)
657 /* i = virtual site | ,k
658 * j = 1st bonded heavy atom | i-j-m
659 * k,l,m = 2nd bonded atoms | `l
663 real bij,bjk,bjl,bjm,aijk,aijl,aijm;
666 bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
667 bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
668 bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
669 bjm = get_bond_length(nrbond, bonds, param->AJ, param->AM);
670 aijk= get_angle (nrang, angles, param->AI, param->AJ, param->AK);
671 aijl= get_angle (nrang, angles, param->AI, param->AJ, param->AL);
672 aijm= get_angle (nrang, angles, param->AI, param->AJ, param->AM);
674 bError = (bij==NOTSET) || (bjk==NOTSET) || (bjl==NOTSET) || (bjm==NOTSET) ||
675 (aijk==NOTSET) || (aijl==NOTSET) || (aijm==NOTSET);
679 /* Calculate component of bond j-k along the direction i-j */
682 /* Calculate component of bond j-l along the direction i-j */
685 /* Calculate component of bond j-m along the direction i-j */
688 if(fabs(pl)<1000*GMX_REAL_MIN || fabs(pm)<1000*GMX_REAL_MIN)
690 fprintf(stderr,"virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
691 param->AI+1,RAD2DEG*aijk,RAD2DEG*aijl,RAD2DEG*aijm);
692 gmx_fatal(FARGS,"invalid construction in calc_vsite4fdn for atom %d: "
693 "pl=%f, pm=%f\n",param->AI+1,pl,pm);
704 fprintf(debug,"params for vsite4fdn %u: %g %g %g\n",
705 param->AI+1,param->C0,param->C1,param->C2);
713 int set_vsites(gmx_bool bVerbose, t_atoms *atoms, gpp_atomtype_t atype,
717 int nvsite,nrbond,nrang,nridih,nrset;
718 gmx_bool bFirst,bSet,bERROR;
719 at2vsitebond_t *at2vb;
728 fprintf(debug, "\nCalculating parameters for virtual sites\n");
730 /* Make a reverse list to avoid ninteractions^2 operations */
731 at2vb = make_at2vsitebond(atoms->nr,plist);
733 for(ftype=0; (ftype<F_NRE); ftype++)
734 if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN) {
736 nvsite+=plist[ftype].nr;
737 for(i=0; (i<plist[ftype].nr); i++) {
738 /* check if all parameters are set */
740 for(j=0; j<NRFP(ftype) && bSet; j++)
741 bSet = plist[ftype].param[i].c[j]!=NOTSET;
744 fprintf(debug,"bSet=%s ",bool_names[bSet]);
745 print_param(debug,ftype,i,&plist[ftype].param[i]);
748 if (bVerbose && bFirst) {
749 fprintf(stderr,"Calculating parameters for virtual sites\n");
753 nrbond=nrang=nridih=0;
758 /* now set the vsite parameters: */
759 get_bondeds(NRAL(ftype), plist[ftype].param[i].a, at2vb, plist,
760 &nrbond, &bonds, &nrang, &angles, &nridih, &idihs);
762 fprintf(debug, "Found %d bonds, %d angles and %d idihs "
763 "for virtual site %u (%s)\n",nrbond,nrang,nridih,
764 plist[ftype].param[i].AI+1,
765 interaction_function[ftype].longname);
766 print_bad(debug, nrbond, bonds, nrang, angles, nridih, idihs);
771 calc_vsite3_param(atype, &(plist[ftype].param[i]), atoms,
772 nrbond, bonds, nrang, angles);
776 calc_vsite3fd_param(&(plist[ftype].param[i]),
777 nrbond, bonds, nrang, angles);
781 calc_vsite3fad_param(&(plist[ftype].param[i]),
782 nrbond, bonds, nrang, angles);
786 calc_vsite3out_param(atype, &(plist[ftype].param[i]), atoms,
787 nrbond, bonds, nrang, angles);
791 calc_vsite4fd_param(&(plist[ftype].param[i]),
792 nrbond, bonds, nrang, angles);
796 calc_vsite4fdn_param(&(plist[ftype].param[i]),
797 nrbond, bonds, nrang, angles);
800 gmx_fatal(FARGS,"Automatic parameter generation not supported "
802 interaction_function[ftype].longname,
803 plist[ftype].param[i].AI+1);
806 gmx_fatal(FARGS,"Automatic parameter generation not supported "
807 "for %s atom %d for this bonding configuration",
808 interaction_function[ftype].longname,
809 plist[ftype].param[i].AI+1);
815 if (debug && plist[ftype].nr)
816 fprintf(stderr,"Calculated parameters for %d out of %d %s atoms\n",
817 nrset,plist[ftype].nr,interaction_function[ftype].longname);
820 done_at2vsitebond(atoms->nr,at2vb);
825 void set_vsites_ptype(gmx_bool bVerbose, gmx_moltype_t *molt)
833 fprintf(stderr,"Setting particle type to V for virtual sites\n");
835 fprintf(stderr,"checking %d functypes\n",F_NRE);
836 for(ftype=0; ftype<F_NRE; ftype++) {
837 il = &molt->ilist[ftype];
838 if (interaction_function[ftype].flags & IF_VSITE) {
839 nra = interaction_function[ftype].nratoms;
844 fprintf(stderr,"doing %d %s virtual sites\n",
845 (nrd / (nra+1)),interaction_function[ftype].longname);
847 for(i=0; (i<nrd); ) {
848 /* The virtual site */
850 molt->atoms.atom[avsite].ptype = eptVSite;
864 static void check_vsite_constraints(t_params *plist,
865 int cftype, int vsite_type[])
872 ps = &(plist[cftype]);
873 for(i=0; (i<ps->nr); i++)
875 atom = ps->param[i].a[k];
876 if (vsite_type[atom]!=NOTSET) {
877 fprintf(stderr,"ERROR: Cannot have constraint (%u-%u) with virtual site (%u)\n",
878 ps->param[i].AI+1, ps->param[i].AJ+1, atom+1);
883 gmx_fatal(FARGS,"There were %d virtual sites involved in constraints",n);
886 static void clean_vsite_bonds(t_params *plist, t_pindex pindex[],
887 int cftype, int vsite_type[])
889 int ftype,i,j,parnr,k,l,m,n,nvsite,nOut,kept_i,vsnral,vsitetype;
890 int nconverted,nremoved;
891 atom_id atom,oatom,constr,at1,at2;
892 atom_id vsiteatoms[MAXATOMLIST];
893 gmx_bool bKeep,bRemove,bUsed,bPresent,bThisFD,bThisOUT,bAllFD,bFirstTwo;
896 if (cftype == F_CONNBONDS)
899 ps = &(plist[cftype]);
905 for(i=0; (i<ps->nr); i++) { /* for all bonds in the plist */
909 /* check if all virtual sites are constructed from the same atoms */
912 fprintf(debug,"constr %u %u:",ps->param[i].AI+1,ps->param[i].AJ+1);
913 for(k=0; (k<2) && !bKeep && !bRemove; k++) {
914 /* for all atoms in the bond */
915 atom = ps->param[i].a[k];
916 if (vsite_type[atom]!=NOTSET) {
918 fprintf(debug," d%d[%d: %d %d %d]",k,atom+1,
919 plist[pindex[atom].ftype].param[pindex[atom].parnr].AJ+1,
920 plist[pindex[atom].ftype].param[pindex[atom].parnr].AK+1,
921 plist[pindex[atom].ftype].param[pindex[atom].parnr].AL+1);
924 bThisFD = ( (pindex[atom].ftype == F_VSITE3FD ) ||
925 (pindex[atom].ftype == F_VSITE3FAD) ||
926 (pindex[atom].ftype == F_VSITE4FD ) ||
927 (pindex[atom].ftype == F_VSITE4FDN ) );
928 bThisOUT= ( (pindex[atom].ftype == F_VSITE3OUT) &&
929 (interaction_function[cftype].flags & IF_CONSTRAINT) );
930 bAllFD = bAllFD && bThisFD;
931 if (bThisFD || bThisOUT) {
932 if(debug)fprintf(debug," %s",bThisOUT?"out":"fd");
933 oatom = ps->param[i].a[1-k]; /* the other atom */
934 if ( vsite_type[oatom]==NOTSET &&
935 oatom==plist[pindex[atom].ftype].param[pindex[atom].parnr].AJ ){
936 /* if the other atom isn't a vsite, and it is AI */
940 if(debug)fprintf(debug," D-AI");
945 /* if this is the first vsite we encounter then
946 store construction atoms */
947 vsnral=NRAL(pindex[atom].ftype)-1;
948 for(m=0; (m<vsnral); m++)
950 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
952 /* if it is not the first then
953 check if this vsite is constructed from the same atoms */
954 if (vsnral == NRAL(pindex[atom].ftype)-1 )
955 for(m=0; (m<vsnral) && !bKeep; m++) {
958 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
959 for(n=0; (n<vsnral) && !bPresent; n++)
960 if (constr == vsiteatoms[n])
964 if(debug)fprintf(debug," !present");
969 if(debug)fprintf(debug," !same#at");
979 /* if we have no virtual sites in this bond, keep it */
981 if (debug)fprintf(debug," no vsite");
985 /* check if all non-vsite atoms are used in construction: */
987 for(k=0; (k<2) && !bKeep; k++) { /* for all atoms in the bond */
988 atom = ps->param[i].a[k];
989 if (vsite_type[atom]==NOTSET) {
991 for(m=0; (m<vsnral) && !bUsed; m++)
992 if (atom == vsiteatoms[m]) {
994 bFirstTwo = bFirstTwo && m<2;
998 if(debug)fprintf(debug," !used");
1003 if ( ! ( bAllFD && bFirstTwo ) )
1004 /* check if all constructing atoms are constrained together */
1005 for (m=0; m<vsnral && !bKeep; m++) { /* all constr. atoms */
1006 at1 = vsiteatoms[m];
1007 at2 = vsiteatoms[(m+1) % vsnral];
1009 for (ftype=0; ftype<F_NRE; ftype++)
1010 if ( interaction_function[ftype].flags & IF_CONSTRAINT )
1011 for (j=0; (j<plist[ftype].nr) && !bPresent; j++)
1012 /* all constraints until one matches */
1013 bPresent = ( ( (plist[ftype].param[j].AI == at1) &&
1014 (plist[ftype].param[j].AJ == at2) ) ||
1015 ( (plist[ftype].param[j].AI == at2) &&
1016 (plist[ftype].param[j].AJ == at1) ) );
1019 if(debug)fprintf(debug," !bonded");
1025 if(debug)fprintf(debug," keeping");
1026 /* now copy the bond to the new array */
1027 ps->param[kept_i] = ps->param[i];
1029 } else if (IS_CHEMBOND(cftype)) {
1030 srenew(plist[F_CONNBONDS].param,plist[F_CONNBONDS].nr+1);
1031 plist[F_CONNBONDS].param[plist[F_CONNBONDS].nr] = ps->param[i];
1032 plist[F_CONNBONDS].nr++;
1036 if(debug)fprintf(debug,"\n");
1040 fprintf(stderr,"Removed %4d %15ss with virtual sites, %5d left\n",
1041 nremoved, interaction_function[cftype].longname, kept_i);
1043 fprintf(stderr,"Converted %4d %15ss with virtual sites to connections, %5d left\n",
1044 nconverted, interaction_function[cftype].longname, kept_i);
1046 fprintf(stderr,"Warning: removed %d %ss with vsite with %s construction\n"
1047 " This vsite construction does not guarantee constant "
1049 " If the constructions were generated by pdb2gmx ignore "
1051 nOut, interaction_function[cftype].longname,
1052 interaction_function[F_VSITE3OUT].longname );
1056 static void clean_vsite_angles(t_params *plist, t_pindex pindex[],
1057 int cftype, int vsite_type[],
1058 at2vsitecon_t *at2vc)
1060 int i,j,parnr,k,l,m,n,nvsite,kept_i,vsnral,vsitetype;
1061 atom_id atom,constr,at1,at2;
1062 atom_id vsiteatoms[MAXATOMLIST];
1063 gmx_bool bKeep,bUsed,bPresent,bAll3FAD,bFirstTwo;
1066 ps = &(plist[cftype]);
1069 for(i=0; (i<ps->nr); i++) { /* for all angles in the plist */
1072 /* check if all virtual sites are constructed from the same atoms */
1074 for(k=0; (k<3) && !bKeep; k++) { /* for all atoms in the angle */
1075 atom = ps->param[i].a[k];
1076 if (vsite_type[atom]!=NOTSET) {
1078 bAll3FAD = bAll3FAD && (pindex[atom].ftype == F_VSITE3FAD);
1080 /* store construction atoms of first vsite */
1081 vsnral=NRAL(pindex[atom].ftype)-1;
1082 for(m=0; (m<vsnral); m++)
1084 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
1086 /* check if this vsite is constructed from the same atoms */
1087 if (vsnral == NRAL(pindex[atom].ftype)-1 )
1088 for(m=0; (m<vsnral) && !bKeep; m++) {
1091 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
1092 for(n=0; (n<vsnral) && !bPresent; n++)
1093 if (constr == vsiteatoms[n])
1103 /* keep all angles with no virtual sites in them or
1104 with virtual sites with more than 3 constr. atoms */
1105 if ( nvsite == 0 && vsnral > 3 )
1108 /* check if all non-vsite atoms are used in construction: */
1110 for(k=0; (k<3) && !bKeep; k++) { /* for all atoms in the angle */
1111 atom = ps->param[i].a[k];
1112 if (vsite_type[atom]==NOTSET) {
1114 for(m=0; (m<vsnral) && !bUsed; m++)
1115 if (atom == vsiteatoms[m]) {
1117 bFirstTwo = bFirstTwo && m<2;
1124 if ( ! ( bAll3FAD && bFirstTwo ) )
1125 /* check if all constructing atoms are constrained together */
1126 for (m=0; m<vsnral && !bKeep; m++) { /* all constr. atoms */
1127 at1 = vsiteatoms[m];
1128 at2 = vsiteatoms[(m+1) % vsnral];
1130 for(j=0; j<at2vc[at1].nr; j++) {
1131 if (at2vc[at1].aj[j] == at2)
1139 /* now copy the angle to the new array */
1140 ps->param[kept_i] = ps->param[i];
1145 if (ps->nr != kept_i)
1146 fprintf(stderr,"Removed %4d %15ss with virtual sites, %5d left\n",
1147 ps->nr-kept_i, interaction_function[cftype].longname, kept_i);
1151 static void clean_vsite_dihs(t_params *plist, t_pindex pindex[],
1152 int cftype, int vsite_type[])
1154 int ftype,i,parnr,k,l,m,n,nvsite,kept_i,vsnral;
1155 atom_id atom,constr;
1156 atom_id vsiteatoms[4];
1157 gmx_bool bKeep,bUsed,bPresent;
1160 ps = &(plist[cftype]);
1164 for(i=0; (i<ps->nr); i++) { /* for all dihedrals in the plist */
1166 /* check if all virtual sites are constructed from the same atoms */
1168 for(k=0; (k<4) && !bKeep; k++) { /* for all atoms in the dihedral */
1169 atom = ps->param[i].a[k];
1170 if (vsite_type[atom]!=NOTSET) {
1173 /* store construction atoms of first vsite */
1174 vsnral=NRAL(pindex[atom].ftype)-1;
1176 for(m=0; (m<vsnral); m++)
1178 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
1180 fprintf(debug,"dih w. vsite: %u %u %u %u\n",
1181 ps->param[i].AI+1,ps->param[i].AJ+1,
1182 ps->param[i].AK+1,ps->param[i].AL+1);
1183 fprintf(debug,"vsite %u from: %u %u %u\n",
1184 atom+1,vsiteatoms[0]+1,vsiteatoms[1]+1,vsiteatoms[2]+1);
1187 /* check if this vsite is constructed from the same atoms */
1188 if (vsnral == NRAL(pindex[atom].ftype)-1 )
1189 for(m=0; (m<vsnral) && !bKeep; m++) {
1192 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
1193 for(n=0; (n<vsnral) && !bPresent; n++)
1194 if (constr == vsiteatoms[n])
1202 /* keep all dihedrals with no virtual sites in them */
1206 /* check if all atoms in dihedral are either virtual sites, or used in
1207 construction of virtual sites. If so, keep it, if not throw away: */
1208 for(k=0; (k<4) && !bKeep; k++) { /* for all atoms in the dihedral */
1209 atom = ps->param[i].a[k];
1210 if (vsite_type[atom]==NOTSET) {
1212 for(m=0; (m<vsnral) && !bUsed; m++)
1213 if (atom == vsiteatoms[m])
1217 if (debug) fprintf(debug,"unused atom in dih: %u\n",atom+1);
1223 ps->param[kept_i] = ps->param[i];
1228 if (ps->nr != kept_i)
1229 fprintf(stderr,"Removed %4d %15ss with virtual sites, %5d left\n",
1230 ps->nr-kept_i, interaction_function[cftype].longname, kept_i);
1234 void clean_vsite_bondeds(t_params *plist, int natoms, gmx_bool bRmVSiteBds)
1236 int i,k,nvsite,ftype,vsite,parnr;
1239 at2vsitecon_t *at2vc;
1241 pindex=0; /* avoid warnings */
1242 /* make vsite_type array */
1243 snew(vsite_type,natoms);
1244 for(i=0; i<natoms; i++)
1245 vsite_type[i]=NOTSET;
1247 for(ftype=0; ftype<F_NRE; ftype++)
1248 if (interaction_function[ftype].flags & IF_VSITE) {
1249 nvsite+=plist[ftype].nr;
1251 while (i < plist[ftype].nr) {
1252 vsite = plist[ftype].param[i].AI;
1253 if ( vsite_type[vsite] == NOTSET)
1254 vsite_type[vsite] = ftype;
1256 gmx_fatal(FARGS,"multiple vsite constructions for atom %d",vsite+1);
1257 if (ftype == F_VSITEN) {
1258 while (i < plist[ftype].nr && plist[ftype].param[i].AI == vsite)
1266 /* the rest only if we have virtual sites: */
1268 fprintf(stderr,"Cleaning up constraints %swith virtual sites\n",
1269 bRmVSiteBds?"and constant bonded interactions ":"");
1271 /* Make a reverse list to avoid ninteractions^2 operations */
1272 at2vc = make_at2vsitecon(natoms,plist);
1274 snew(pindex,natoms);
1275 for(ftype=0; ftype<F_NRE; ftype++) {
1276 if ((interaction_function[ftype].flags & IF_VSITE) &&
1277 ftype != F_VSITEN) {
1278 for (parnr=0; (parnr<plist[ftype].nr); parnr++) {
1279 k=plist[ftype].param[parnr].AI;
1280 pindex[k].ftype=ftype;
1281 pindex[k].parnr=parnr;
1287 for(i=0; i<natoms; i++)
1288 fprintf(debug,"atom %d vsite_type %s\n",i,
1289 vsite_type[i]==NOTSET ? "NOTSET" :
1290 interaction_function[vsite_type[i]].name);
1292 /* remove things with vsite atoms */
1293 for(ftype=0; ftype<F_NRE; ftype++)
1294 if ( ( ( interaction_function[ftype].flags & IF_BOND ) && bRmVSiteBds ) ||
1295 ( interaction_function[ftype].flags & IF_CONSTRAINT ) ) {
1296 if (interaction_function[ftype].flags & (IF_BTYPE | IF_CONSTRAINT) )
1297 clean_vsite_bonds (plist, pindex, ftype, vsite_type);
1298 else if (interaction_function[ftype].flags & IF_ATYPE)
1299 clean_vsite_angles(plist, pindex, ftype, vsite_type, at2vc);
1300 else if ( (ftype==F_PDIHS) || (ftype==F_IDIHS) )
1301 clean_vsite_dihs (plist, pindex, ftype, vsite_type);
1303 /* check if we have constraints left with virtual sites in them */
1304 for(ftype=0; ftype<F_NRE; ftype++)
1305 if (interaction_function[ftype].flags & IF_CONSTRAINT)
1306 check_vsite_constraints(plist, ftype, vsite_type);
1308 done_at2vsitecon(natoms,at2vc);