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
42 #include "vsite_parm.h"
51 #include "gmx_fatal.h"
68 vsitebondparam_t *vsbp;
76 static int vsite_bond_nrcheck(int ftype)
80 if ((interaction_function[ftype].flags & (IF_BTYPE | IF_CONSTRAINT | IF_ATYPE)) || (ftype == F_IDIHS))
81 nrcheck = NRAL(ftype);
88 static void enter_bonded(int nratoms, int *nrbonded, t_mybonded **bondeds,
93 srenew(*bondeds, *nrbonded+1);
95 /* copy atom numbers */
96 for(j=0; j<nratoms; j++)
97 (*bondeds)[*nrbonded].a[j] = param->a[j];
99 (*bondeds)[*nrbonded].c = param->C0;
104 static void get_bondeds(int nrat, t_iatom atoms[],
105 at2vsitebond_t *at2vb, t_params plist[],
106 int *nrbond, t_mybonded **bonds,
107 int *nrang, t_mybonded **angles,
108 int *nridih, t_mybonded **idihs )
110 int k,i,ftype,nrcheck;
113 for(k=0; k<nrat; k++) {
114 for(i=0; i<at2vb[atoms[k]].nr; i++) {
115 ftype = at2vb[atoms[k]].vsbp[i].ftype;
116 param = at2vb[atoms[k]].vsbp[i].param;
117 nrcheck = vsite_bond_nrcheck(ftype);
118 /* abuse nrcheck to see if we're adding bond, angle or idih */
120 case 2: enter_bonded(nrcheck,nrbond,bonds, param); break;
121 case 3: enter_bonded(nrcheck,nrang, angles,param); break;
122 case 4: enter_bonded(nrcheck,nridih,idihs, param); break;
128 static at2vsitebond_t *make_at2vsitebond(int natoms,t_params plist[])
131 int ftype,i,j,nrcheck,nr;
133 at2vsitebond_t *at2vb;
138 for(ftype=0; (ftype<F_NRE); ftype++) {
139 if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN) {
140 for(i=0; (i<plist[ftype].nr); i++) {
141 for(j=0; j<NRAL(ftype); j++)
142 bVSI[plist[ftype].param[i].a[j]] = TRUE;
147 for(ftype=0; (ftype<F_NRE); ftype++) {
148 nrcheck = vsite_bond_nrcheck(ftype);
150 for(i=0; (i<plist[ftype].nr); i++) {
151 aa = plist[ftype].param[i].a;
152 for(j=0; j<nrcheck; j++) {
154 nr = at2vb[aa[j]].nr;
156 srenew(at2vb[aa[j]].vsbp,nr+10);
157 at2vb[aa[j]].vsbp[nr].ftype = ftype;
158 at2vb[aa[j]].vsbp[nr].param = &plist[ftype].param[i];
170 static void done_at2vsitebond(int natoms,at2vsitebond_t *at2vb)
174 for(i=0; i<natoms; i++)
176 sfree(at2vb[i].vsbp);
180 static at2vsitecon_t *make_at2vsitecon(int natoms,t_params plist[])
183 int ftype,i,j,ai,aj,nr;
184 at2vsitecon_t *at2vc;
189 for(ftype=0; (ftype<F_NRE); ftype++) {
190 if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN) {
191 for(i=0; (i<plist[ftype].nr); i++) {
192 for(j=0; j<NRAL(ftype); j++)
193 bVSI[plist[ftype].param[i].a[j]] = TRUE;
198 for(ftype=0; (ftype<F_NRE); ftype++) {
199 if (interaction_function[ftype].flags & IF_CONSTRAINT) {
200 for(i=0; (i<plist[ftype].nr); i++) {
201 ai = plist[ftype].param[i].AI;
202 aj = plist[ftype].param[i].AJ;
203 if (bVSI[ai] && bVSI[aj]) {
204 /* Store forward direction */
207 srenew(at2vc[ai].aj,nr+10);
208 at2vc[ai].aj[nr] = aj;
210 /* Store backward direction */
213 srenew(at2vc[aj].aj,nr+10);
214 at2vc[aj].aj[nr] = ai;
225 static void done_at2vsitecon(int natoms,at2vsitecon_t *at2vc)
229 for(i=0; i<natoms; i++)
236 static void print_bad(FILE *fp,
237 int nrbond, t_mybonded *bonds,
238 int nrang, t_mybonded *angles,
239 int nridih, t_mybonded *idihs )
244 fprintf(fp,"bonds:");
245 for(i=0; i<nrbond; i++)
246 fprintf(fp," %u-%u (%g)",
247 bonds[i].AI+1, bonds[i].AJ+1, bonds[i].c);
251 fprintf(fp,"angles:");
252 for(i=0; i<nrang; i++)
253 fprintf(fp," %u-%u-%u (%g)",
254 angles[i].AI+1, angles[i].AJ+1,
255 angles[i].AK+1, angles[i].c);
259 fprintf(fp,"idihs:");
260 for(i=0; i<nridih; i++)
261 fprintf(fp," %u-%u-%u-%u (%g)",
262 idihs[i].AI+1, idihs[i].AJ+1,
263 idihs[i].AK+1, idihs[i].AL+1, idihs[i].c);
268 static void print_param(FILE *fp, int ftype, int i, t_param *param)
271 static int prev_ftype= NOTSET;
272 static int prev_i = NOTSET;
275 if ( (ftype!=prev_ftype) || (i!=prev_i) ) {
280 fprintf(fp,"(%d) plist[%s].param[%d]",
281 pass,interaction_function[ftype].name,i);
282 for(j=0; j<NRFP(ftype); j++)
283 fprintf(fp,".c[%d]=%g ",j,param->c[j]);
288 static real get_bond_length(int nrbond, t_mybonded bonds[],
289 t_iatom ai, t_iatom aj)
295 for (i=0; i < nrbond && (bondlen==NOTSET); i++) {
296 /* check both ways */
297 if ( ( (ai == bonds[i].AI) && (aj == bonds[i].AJ) ) ||
298 ( (ai == bonds[i].AJ) && (aj == bonds[i].AI) ) )
299 bondlen = bonds[i].c; /* note: bonds[i].c might be NOTSET */
304 static real get_angle(int nrang, t_mybonded angles[],
305 t_iatom ai, t_iatom aj, t_iatom ak)
311 for (i=0; i < nrang && (angle==NOTSET); i++) {
312 /* check both ways */
313 if ( ( (ai==angles[i].AI) && (aj==angles[i].AJ) && (ak==angles[i].AK) ) ||
314 ( (ai==angles[i].AK) && (aj==angles[i].AJ) && (ak==angles[i].AI) ) )
315 angle = DEG2RAD*angles[i].c;
320 static char *get_atomtype_name_AB(t_atom *atom,gpp_atomtype_t atype)
324 name = get_atomtype_name(atom->type,atype);
326 /* When using the decoupling option, atom types are changed
327 * to decoupled for the non-bonded interactions, but the virtual
328 * sites constructions should be based on the "normal" interactions.
329 * So we return the state B atom type name if the state A atom
330 * type is the decoupled one. We should actually check for the atom
331 * type number, but that's not passed here. So we check for
332 * the decoupled atom type name. This should not cause trouble
333 * as this code is only used for topologies with v-sites without
334 * parameters generated by pdb2gmx.
336 if (strcmp(name,"decoupled") == 0)
338 name = get_atomtype_name(atom->typeB,atype);
344 static gmx_bool calc_vsite3_param(gpp_atomtype_t atype,
345 t_param *param, t_atoms *at,
346 int nrbond, t_mybonded *bonds,
347 int nrang, t_mybonded *angles )
349 /* i = virtual site | ,k
350 * j = 1st bonded heavy atom | i-j
351 * k,l = 2nd bonded atoms | `l
354 gmx_bool bXH3,bError;
355 real bjk,bjl,a=-1,b=-1;
356 /* check if this is part of a NH3 , NH2-umbrella or CH3 group,
357 * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
361 fprintf(debug,"atom %u type %s ",
363 get_atomtype_name_AB(&at->atom[param->a[i]],atype));
367 ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AK],atype),"MNH",3)==0) &&
368 (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AL],atype),"MNH",3)==0) ) ||
369 ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AK],atype),"MCH3",4)==0) &&
370 (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AL],atype),"MCH3",4)==0) );
372 bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
373 bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
374 bError = (bjk==NOTSET) || (bjl==NOTSET);
376 /* now we get some XH2/XH3 group specific construction */
377 /* note: we call the heavy atom 'C' and the X atom 'N' */
378 real bMM,bCM,bCN,bNH,aCNH,dH,rH,dM,rM;
381 /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
382 bError = bError || (bjk!=bjl);
384 /* the X atom (C or N) in the XH2/XH3 group is the first after the masses: */
385 aN = max(param->AK,param->AL)+1;
387 /* get common bonds */
388 bMM = get_bond_length(nrbond, bonds, param->AK, param->AL);
390 bCN = get_bond_length(nrbond, bonds, param->AJ, aN);
391 bError = bError || (bMM==NOTSET) || (bCN==NOTSET);
393 /* calculate common things */
395 dM = sqrt( sqr(bCM) - sqr(rM) );
397 /* are we dealing with the X atom? */
398 if ( param->AI == aN ) {
399 /* this is trivial */
400 a = b = 0.5 * bCN/dM;
403 /* get other bondlengths and angles: */
404 bNH = get_bond_length(nrbond, bonds, aN, param->AI);
405 aCNH= get_angle (nrang, angles, param->AJ, aN, param->AI);
406 bError = bError || (bNH==NOTSET) || (aCNH==NOTSET);
409 dH = bCN - bNH * cos(aCNH);
410 rH = bNH * sin(aCNH);
412 a = 0.5 * ( dH/dM + rH/rM );
413 b = 0.5 * ( dH/dM - rH/rM );
416 gmx_fatal(FARGS,"calc_vsite3_param not implemented for the general case "
417 "(atom %d)",param->AI+1);
423 fprintf(debug,"params for vsite3 %u: %g %g\n",
424 param->AI+1,param->C0,param->C1);
429 static gmx_bool calc_vsite3fd_param(t_param *param,
430 int nrbond, t_mybonded *bonds,
431 int nrang, t_mybonded *angles)
433 /* i = virtual site | ,k
434 * j = 1st bonded heavy atom | i-j
435 * k,l = 2nd bonded atoms | `l
439 real bij,bjk,bjl,aijk,aijl,rk,rl;
441 bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
442 bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
443 bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
444 aijk= get_angle (nrang, angles, param->AI, param->AJ, param->AK);
445 aijl= get_angle (nrang, angles, param->AI, param->AJ, param->AL);
446 bError = (bij==NOTSET) || (bjk==NOTSET) || (bjl==NOTSET) ||
447 (aijk==NOTSET) || (aijl==NOTSET);
449 rk = bjk * sin(aijk);
450 rl = bjl * sin(aijl);
451 param->C0 = rk / (rk + rl);
452 param->C1 = -bij; /* 'bond'-length for fixed distance vsite */
455 fprintf(debug,"params for vsite3fd %u: %g %g\n",
456 param->AI+1,param->C0,param->C1);
460 static gmx_bool calc_vsite3fad_param(t_param *param,
461 int nrbond, t_mybonded *bonds,
462 int nrang, t_mybonded *angles)
464 /* i = virtual site |
465 * j = 1st bonded heavy atom | i-j
466 * k = 2nd bonded heavy atom | `k-l
467 * l = 3d bonded heavy atom |
470 gmx_bool bSwapParity,bError;
473 bSwapParity = ( param->C1 == -1 );
475 bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
476 aijk = get_angle (nrang, angles, param->AI, param->AJ, param->AK);
477 bError = (bij==NOTSET) || (aijk==NOTSET);
479 param->C1 = bij; /* 'bond'-length for fixed distance vsite */
480 param->C0 = RAD2DEG*aijk; /* 'bond'-angle for fixed angle vsite */
483 param->C0 = 360 - param->C0;
486 fprintf(debug,"params for vsite3fad %u: %g %g\n",
487 param->AI+1,param->C0,param->C1);
491 static gmx_bool calc_vsite3out_param(gpp_atomtype_t atype,
492 t_param *param, t_atoms *at,
493 int nrbond, t_mybonded *bonds,
494 int nrang, t_mybonded *angles)
496 /* i = virtual site | ,k
497 * j = 1st bonded heavy atom | i-j
498 * k,l = 2nd bonded atoms | `l
499 * NOTE: i is out of the j-k-l plane!
502 gmx_bool bXH3,bError,bSwapParity;
503 real bij,bjk,bjl,aijk,aijl,akjl,pijk,pijl,a,b,c;
505 /* check if this is part of a NH2-umbrella, NH3 or CH3 group,
506 * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
510 fprintf(debug,"atom %u type %s ",
511 param->a[i]+1,get_atomtype_name_AB(&at->atom[param->a[i]],atype));
515 ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AK],atype),"MNH",3)==0) &&
516 (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AL],atype),"MNH",3)==0) ) ||
517 ( (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AK],atype),"MCH3",4)==0) &&
518 (gmx_strncasecmp(get_atomtype_name_AB(&at->atom[param->AL],atype),"MCH3",4)==0) );
520 /* check if construction parity must be swapped */
521 bSwapParity = ( param->C1 == -1 );
523 bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
524 bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
525 bError = (bjk==NOTSET) || (bjl==NOTSET);
527 /* now we get some XH3 group specific construction */
528 /* note: we call the heavy atom 'C' and the X atom 'N' */
529 real bMM,bCM,bCN,bNH,aCNH,dH,rH,rHx,rHy,dM,rM;
532 /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
533 bError = bError || (bjk!=bjl);
535 /* the X atom (C or N) in the XH3 group is the first after the masses: */
536 aN = max(param->AK,param->AL)+1;
538 /* get all bondlengths and angles: */
539 bMM = get_bond_length(nrbond, bonds, param->AK, param->AL);
541 bCN = get_bond_length(nrbond, bonds, param->AJ, aN);
542 bNH = get_bond_length(nrbond, bonds, aN, param->AI);
543 aCNH= get_angle (nrang, angles, param->AJ, aN, param->AI);
545 (bMM==NOTSET) || (bCN==NOTSET) || (bNH==NOTSET) || (aCNH==NOTSET);
548 dH = bCN - bNH * cos(aCNH);
549 rH = bNH * sin(aCNH);
550 /* we assume the H's are symmetrically distributed */
551 rHx = rH*cos(DEG2RAD*30);
552 rHy = rH*sin(DEG2RAD*30);
554 dM = sqrt( sqr(bCM) - sqr(rM) );
555 a = 0.5*( (dH/dM) - (rHy/rM) );
556 b = 0.5*( (dH/dM) + (rHy/rM) );
560 /* this is the general construction */
562 bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
563 aijk= get_angle (nrang, angles, param->AI, param->AJ, param->AK);
564 aijl= get_angle (nrang, angles, param->AI, param->AJ, param->AL);
565 akjl= get_angle (nrang, angles, param->AK, param->AJ, param->AL);
567 (bij==NOTSET) || (aijk==NOTSET) || (aijl==NOTSET) || (akjl==NOTSET);
569 pijk = cos(aijk)*bij;
570 pijl = cos(aijl)*bij;
571 a = ( pijk + (pijk*cos(akjl)-pijl) * cos(akjl) / sqr(sin(akjl)) ) / bjk;
572 b = ( pijl + (pijl*cos(akjl)-pijk) * cos(akjl) / sqr(sin(akjl)) ) / bjl;
573 c = - sqrt( sqr(bij) -
574 ( sqr(pijk) - 2*pijk*pijl*cos(akjl) + sqr(pijl) )
576 / ( bjk*bjl*sin(akjl) );
586 fprintf(debug,"params for vsite3out %u: %g %g %g\n",
587 param->AI+1,param->C0,param->C1,param->C2);
591 static gmx_bool calc_vsite4fd_param(t_param *param,
592 int nrbond, t_mybonded *bonds,
593 int nrang, t_mybonded *angles)
595 /* i = virtual site | ,k
596 * j = 1st bonded heavy atom | i-j-m
597 * k,l,m = 2nd bonded atoms | `l
601 real bij,bjk,bjl,bjm,aijk,aijl,aijm,akjm,akjl;
602 real pk,pl,pm,cosakl,cosakm,sinakl,sinakm,cl,cm;
604 bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
605 bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
606 bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
607 bjm = get_bond_length(nrbond, bonds, param->AJ, param->AM);
608 aijk= get_angle (nrang, angles, param->AI, param->AJ, param->AK);
609 aijl= get_angle (nrang, angles, param->AI, param->AJ, param->AL);
610 aijm= get_angle (nrang, angles, param->AI, param->AJ, param->AM);
611 akjm= get_angle (nrang, angles, param->AK, param->AJ, param->AM);
612 akjl= get_angle (nrang, angles, param->AK, param->AJ, param->AL);
613 bError = (bij==NOTSET) || (bjk==NOTSET) || (bjl==NOTSET) || (bjm==NOTSET) ||
614 (aijk==NOTSET) || (aijl==NOTSET) || (aijm==NOTSET) || (akjm==NOTSET) ||
621 cosakl = (cos(akjl) - cos(aijk)*cos(aijl)) / (sin(aijk)*sin(aijl));
622 cosakm = (cos(akjm) - cos(aijk)*cos(aijm)) / (sin(aijk)*sin(aijm));
623 if ( cosakl < -1 || cosakl > 1 || cosakm < -1 || cosakm > 1 ) {
624 fprintf(stderr,"virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
625 param->AI+1,RAD2DEG*aijk,RAD2DEG*aijl,RAD2DEG*aijm);
626 gmx_fatal(FARGS,"invalid construction in calc_vsite4fd for atom %d: "
627 "cosakl=%f, cosakm=%f\n",param->AI+1,cosakl,cosakm);
629 sinakl = sqrt(1-sqr(cosakl));
630 sinakm = sqrt(1-sqr(cosakm));
632 /* note: there is a '+' because of the way the sines are calculated */
633 cl = -pk / ( pl*cosakl - pk + pl*sinakl*(pm*cosakm-pk)/(pm*sinakm) );
634 cm = -pk / ( pm*cosakm - pk + pm*sinakm*(pl*cosakl-pk)/(pl*sinakl) );
640 fprintf(debug,"params for vsite4fd %u: %g %g %g\n",
641 param->AI+1,param->C0,param->C1,param->C2);
649 calc_vsite4fdn_param(t_param *param,
650 int nrbond, t_mybonded *bonds,
651 int nrang, t_mybonded *angles)
653 /* i = virtual site | ,k
654 * j = 1st bonded heavy atom | i-j-m
655 * k,l,m = 2nd bonded atoms | `l
659 real bij,bjk,bjl,bjm,aijk,aijl,aijm;
662 bij = get_bond_length(nrbond, bonds, param->AI, param->AJ);
663 bjk = get_bond_length(nrbond, bonds, param->AJ, param->AK);
664 bjl = get_bond_length(nrbond, bonds, param->AJ, param->AL);
665 bjm = get_bond_length(nrbond, bonds, param->AJ, param->AM);
666 aijk= get_angle (nrang, angles, param->AI, param->AJ, param->AK);
667 aijl= get_angle (nrang, angles, param->AI, param->AJ, param->AL);
668 aijm= get_angle (nrang, angles, param->AI, param->AJ, param->AM);
670 bError = (bij==NOTSET) || (bjk==NOTSET) || (bjl==NOTSET) || (bjm==NOTSET) ||
671 (aijk==NOTSET) || (aijl==NOTSET) || (aijm==NOTSET);
675 /* Calculate component of bond j-k along the direction i-j */
678 /* Calculate component of bond j-l along the direction i-j */
681 /* Calculate component of bond j-m along the direction i-j */
684 if(fabs(pl)<1000*GMX_REAL_MIN || fabs(pm)<1000*GMX_REAL_MIN)
686 fprintf(stderr,"virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
687 param->AI+1,RAD2DEG*aijk,RAD2DEG*aijl,RAD2DEG*aijm);
688 gmx_fatal(FARGS,"invalid construction in calc_vsite4fdn for atom %d: "
689 "pl=%f, pm=%f\n",param->AI+1,pl,pm);
700 fprintf(debug,"params for vsite4fdn %u: %g %g %g\n",
701 param->AI+1,param->C0,param->C1,param->C2);
709 int set_vsites(gmx_bool bVerbose, t_atoms *atoms, gpp_atomtype_t atype,
713 int nvsite,nrbond,nrang,nridih,nrset;
714 gmx_bool bFirst,bSet,bERROR;
715 at2vsitebond_t *at2vb;
724 fprintf(debug, "\nCalculating parameters for virtual sites\n");
726 /* Make a reverse list to avoid ninteractions^2 operations */
727 at2vb = make_at2vsitebond(atoms->nr,plist);
729 for(ftype=0; (ftype<F_NRE); ftype++)
730 if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN) {
732 nvsite+=plist[ftype].nr;
733 for(i=0; (i<plist[ftype].nr); i++) {
734 /* check if all parameters are set */
736 for(j=0; j<NRFP(ftype) && bSet; j++)
737 bSet = plist[ftype].param[i].c[j]!=NOTSET;
740 fprintf(debug,"bSet=%s ",bool_names[bSet]);
741 print_param(debug,ftype,i,&plist[ftype].param[i]);
744 if (bVerbose && bFirst) {
745 fprintf(stderr,"Calculating parameters for virtual sites\n");
749 nrbond=nrang=nridih=0;
754 /* now set the vsite parameters: */
755 get_bondeds(NRAL(ftype), plist[ftype].param[i].a, at2vb, plist,
756 &nrbond, &bonds, &nrang, &angles, &nridih, &idihs);
758 fprintf(debug, "Found %d bonds, %d angles and %d idihs "
759 "for virtual site %u (%s)\n",nrbond,nrang,nridih,
760 plist[ftype].param[i].AI+1,
761 interaction_function[ftype].longname);
762 print_bad(debug, nrbond, bonds, nrang, angles, nridih, idihs);
767 calc_vsite3_param(atype, &(plist[ftype].param[i]), atoms,
768 nrbond, bonds, nrang, angles);
772 calc_vsite3fd_param(&(plist[ftype].param[i]),
773 nrbond, bonds, nrang, angles);
777 calc_vsite3fad_param(&(plist[ftype].param[i]),
778 nrbond, bonds, nrang, angles);
782 calc_vsite3out_param(atype, &(plist[ftype].param[i]), atoms,
783 nrbond, bonds, nrang, angles);
787 calc_vsite4fd_param(&(plist[ftype].param[i]),
788 nrbond, bonds, nrang, angles);
792 calc_vsite4fdn_param(&(plist[ftype].param[i]),
793 nrbond, bonds, nrang, angles);
796 gmx_fatal(FARGS,"Automatic parameter generation not supported "
798 interaction_function[ftype].longname,
799 plist[ftype].param[i].AI+1);
802 gmx_fatal(FARGS,"Automatic parameter generation not supported "
803 "for %s atom %d for this bonding configuration",
804 interaction_function[ftype].longname,
805 plist[ftype].param[i].AI+1);
811 if (debug && plist[ftype].nr)
812 fprintf(stderr,"Calculated parameters for %d out of %d %s atoms\n",
813 nrset,plist[ftype].nr,interaction_function[ftype].longname);
816 done_at2vsitebond(atoms->nr,at2vb);
821 void set_vsites_ptype(gmx_bool bVerbose, gmx_moltype_t *molt)
829 fprintf(stderr,"Setting particle type to V for virtual sites\n");
831 fprintf(stderr,"checking %d functypes\n",F_NRE);
832 for(ftype=0; ftype<F_NRE; ftype++) {
833 il = &molt->ilist[ftype];
834 if (interaction_function[ftype].flags & IF_VSITE) {
835 nra = interaction_function[ftype].nratoms;
840 fprintf(stderr,"doing %d %s virtual sites\n",
841 (nrd / (nra+1)),interaction_function[ftype].longname);
843 for(i=0; (i<nrd); ) {
844 /* The virtual site */
846 molt->atoms.atom[avsite].ptype = eptVSite;
860 static void check_vsite_constraints(t_params *plist,
861 int cftype, int vsite_type[])
868 ps = &(plist[cftype]);
869 for(i=0; (i<ps->nr); i++)
871 atom = ps->param[i].a[k];
872 if (vsite_type[atom]!=NOTSET) {
873 fprintf(stderr,"ERROR: Cannot have constraint (%u-%u) with virtual site (%u)\n",
874 ps->param[i].AI+1, ps->param[i].AJ+1, atom+1);
879 gmx_fatal(FARGS,"There were %d virtual sites involved in constraints",n);
882 static void clean_vsite_bonds(t_params *plist, t_pindex pindex[],
883 int cftype, int vsite_type[])
885 int ftype,i,j,parnr,k,l,m,n,nvsite,nOut,kept_i,vsnral,vsitetype;
886 int nconverted,nremoved;
887 atom_id atom,oatom,constr,at1,at2;
888 atom_id vsiteatoms[MAXATOMLIST];
889 gmx_bool bKeep,bRemove,bUsed,bPresent,bThisFD,bThisOUT,bAllFD,bFirstTwo;
892 if (cftype == F_CONNBONDS)
895 ps = &(plist[cftype]);
901 for(i=0; (i<ps->nr); i++) { /* for all bonds in the plist */
905 /* check if all virtual sites are constructed from the same atoms */
908 fprintf(debug,"constr %u %u:",ps->param[i].AI+1,ps->param[i].AJ+1);
909 for(k=0; (k<2) && !bKeep && !bRemove; k++) {
910 /* for all atoms in the bond */
911 atom = ps->param[i].a[k];
912 if (vsite_type[atom]!=NOTSET) {
914 fprintf(debug," d%d[%d: %d %d %d]",k,atom+1,
915 plist[pindex[atom].ftype].param[pindex[atom].parnr].AJ+1,
916 plist[pindex[atom].ftype].param[pindex[atom].parnr].AK+1,
917 plist[pindex[atom].ftype].param[pindex[atom].parnr].AL+1);
920 bThisFD = ( (pindex[atom].ftype == F_VSITE3FD ) ||
921 (pindex[atom].ftype == F_VSITE3FAD) ||
922 (pindex[atom].ftype == F_VSITE4FD ) ||
923 (pindex[atom].ftype == F_VSITE4FDN ) );
924 bThisOUT= ( (pindex[atom].ftype == F_VSITE3OUT) &&
925 (interaction_function[cftype].flags & IF_CONSTRAINT) );
926 bAllFD = bAllFD && bThisFD;
927 if (bThisFD || bThisOUT) {
928 if(debug)fprintf(debug," %s",bThisOUT?"out":"fd");
929 oatom = ps->param[i].a[1-k]; /* the other atom */
930 if ( vsite_type[oatom]==NOTSET &&
931 oatom==plist[pindex[atom].ftype].param[pindex[atom].parnr].AJ ){
932 /* if the other atom isn't a vsite, and it is AI */
936 if(debug)fprintf(debug," D-AI");
941 /* if this is the first vsite we encounter then
942 store construction atoms */
943 vsnral=NRAL(pindex[atom].ftype)-1;
944 for(m=0; (m<vsnral); m++)
946 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
948 /* if it is not the first then
949 check if this vsite is constructed from the same atoms */
950 if (vsnral == NRAL(pindex[atom].ftype)-1 )
951 for(m=0; (m<vsnral) && !bKeep; m++) {
954 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
955 for(n=0; (n<vsnral) && !bPresent; n++)
956 if (constr == vsiteatoms[n])
960 if(debug)fprintf(debug," !present");
965 if(debug)fprintf(debug," !same#at");
975 /* if we have no virtual sites in this bond, keep it */
977 if (debug)fprintf(debug," no vsite");
981 /* check if all non-vsite atoms are used in construction: */
983 for(k=0; (k<2) && !bKeep; k++) { /* for all atoms in the bond */
984 atom = ps->param[i].a[k];
985 if (vsite_type[atom]==NOTSET) {
987 for(m=0; (m<vsnral) && !bUsed; m++)
988 if (atom == vsiteatoms[m]) {
990 bFirstTwo = bFirstTwo && m<2;
994 if(debug)fprintf(debug," !used");
999 if ( ! ( bAllFD && bFirstTwo ) )
1000 /* check if all constructing atoms are constrained together */
1001 for (m=0; m<vsnral && !bKeep; m++) { /* all constr. atoms */
1002 at1 = vsiteatoms[m];
1003 at2 = vsiteatoms[(m+1) % vsnral];
1005 for (ftype=0; ftype<F_NRE; ftype++)
1006 if ( interaction_function[ftype].flags & IF_CONSTRAINT )
1007 for (j=0; (j<plist[ftype].nr) && !bPresent; j++)
1008 /* all constraints until one matches */
1009 bPresent = ( ( (plist[ftype].param[j].AI == at1) &&
1010 (plist[ftype].param[j].AJ == at2) ) ||
1011 ( (plist[ftype].param[j].AI == at2) &&
1012 (plist[ftype].param[j].AJ == at1) ) );
1015 if(debug)fprintf(debug," !bonded");
1021 if(debug)fprintf(debug," keeping");
1022 /* now copy the bond to the new array */
1023 memcpy(&(ps->param[kept_i]),
1024 &(ps->param[i]),(size_t)sizeof(ps->param[0]));
1026 } else if (IS_CHEMBOND(cftype)) {
1027 srenew(plist[F_CONNBONDS].param,plist[F_CONNBONDS].nr+1);
1028 memcpy(&(plist[F_CONNBONDS].param[plist[F_CONNBONDS].nr]),
1029 &(ps->param[i]),(size_t)sizeof(plist[F_CONNBONDS].param[0]));
1030 plist[F_CONNBONDS].nr++;
1034 if(debug)fprintf(debug,"\n");
1038 fprintf(stderr,"Removed %4d %15ss with virtual sites, %5d left\n",
1039 nremoved, interaction_function[cftype].longname, kept_i);
1041 fprintf(stderr,"Converted %4d %15ss with virtual sites to connections, %5d left\n",
1042 nconverted, interaction_function[cftype].longname, kept_i);
1044 fprintf(stderr,"Warning: removed %d %ss with vsite with %s construction\n"
1045 " This vsite construction does not guarantee constant "
1047 " If the constructions were generated by pdb2gmx ignore "
1049 nOut, interaction_function[cftype].longname,
1050 interaction_function[F_VSITE3OUT].longname );
1054 static void clean_vsite_angles(t_params *plist, t_pindex pindex[],
1055 int cftype, int vsite_type[],
1056 at2vsitecon_t *at2vc)
1058 int i,j,parnr,k,l,m,n,nvsite,kept_i,vsnral,vsitetype;
1059 atom_id atom,constr,at1,at2;
1060 atom_id vsiteatoms[MAXATOMLIST];
1061 gmx_bool bKeep,bUsed,bPresent,bAll3FAD,bFirstTwo;
1064 ps = &(plist[cftype]);
1067 for(i=0; (i<ps->nr); i++) { /* for all angles in the plist */
1070 /* check if all virtual sites are constructed from the same atoms */
1072 for(k=0; (k<3) && !bKeep; k++) { /* for all atoms in the angle */
1073 atom = ps->param[i].a[k];
1074 if (vsite_type[atom]!=NOTSET) {
1076 bAll3FAD = bAll3FAD && (pindex[atom].ftype == F_VSITE3FAD);
1078 /* store construction atoms of first vsite */
1079 vsnral=NRAL(pindex[atom].ftype)-1;
1080 for(m=0; (m<vsnral); m++)
1082 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
1084 /* check if this vsite is constructed from the same atoms */
1085 if (vsnral == NRAL(pindex[atom].ftype)-1 )
1086 for(m=0; (m<vsnral) && !bKeep; m++) {
1089 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
1090 for(n=0; (n<vsnral) && !bPresent; n++)
1091 if (constr == vsiteatoms[n])
1101 /* keep all angles with no virtual sites in them or
1102 with virtual sites with more than 3 constr. atoms */
1103 if ( nvsite == 0 && vsnral > 3 )
1106 /* check if all non-vsite atoms are used in construction: */
1108 for(k=0; (k<3) && !bKeep; k++) { /* for all atoms in the angle */
1109 atom = ps->param[i].a[k];
1110 if (vsite_type[atom]==NOTSET) {
1112 for(m=0; (m<vsnral) && !bUsed; m++)
1113 if (atom == vsiteatoms[m]) {
1115 bFirstTwo = bFirstTwo && m<2;
1122 if ( ! ( bAll3FAD && bFirstTwo ) )
1123 /* check if all constructing atoms are constrained together */
1124 for (m=0; m<vsnral && !bKeep; m++) { /* all constr. atoms */
1125 at1 = vsiteatoms[m];
1126 at2 = vsiteatoms[(m+1) % vsnral];
1128 for(j=0; j<at2vc[at1].nr; j++) {
1129 if (at2vc[at1].aj[j] == at2)
1137 /* now copy the angle to the new array */
1138 memcpy(&(ps->param[kept_i]),
1139 &(ps->param[i]),(size_t)sizeof(ps->param[0]));
1144 if (ps->nr != kept_i)
1145 fprintf(stderr,"Removed %4d %15ss with virtual sites, %5d left\n",
1146 ps->nr-kept_i, interaction_function[cftype].longname, kept_i);
1150 static void clean_vsite_dihs(t_params *plist, t_pindex pindex[],
1151 int cftype, int vsite_type[])
1153 int ftype,i,parnr,k,l,m,n,nvsite,kept_i,vsnral;
1154 atom_id atom,constr;
1155 atom_id vsiteatoms[3];
1156 gmx_bool bKeep,bUsed,bPresent;
1159 ps = &(plist[cftype]);
1163 for(i=0; (i<ps->nr); i++) { /* for all dihedrals in the plist */
1165 /* check if all virtual sites are constructed from the same atoms */
1167 for(k=0; (k<4) && !bKeep; k++) { /* for all atoms in the dihedral */
1168 atom = ps->param[i].a[k];
1169 if (vsite_type[atom]!=NOTSET) {
1172 /* store construction atoms of first vsite */
1173 vsnral=NRAL(pindex[atom].ftype)-1;
1174 for(m=0; (m<vsnral); m++)
1176 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
1178 fprintf(debug,"dih w. vsite: %u %u %u %u\n",
1179 ps->param[i].AI+1,ps->param[i].AJ+1,
1180 ps->param[i].AK+1,ps->param[i].AL+1);
1181 fprintf(debug,"vsite %u from: %u %u %u\n",
1182 atom+1,vsiteatoms[0]+1,vsiteatoms[1]+1,vsiteatoms[2]+1);
1185 /* check if this vsite is constructed from the same atoms */
1186 if (vsnral == NRAL(pindex[atom].ftype)-1 )
1187 for(m=0; (m<vsnral) && !bKeep; m++) {
1190 plist[pindex[atom].ftype].param[pindex[atom].parnr].a[m+1];
1191 for(n=0; (n<vsnral) && !bPresent; n++)
1192 if (constr == vsiteatoms[n])
1200 /* keep all dihedrals with no virtual sites in them */
1204 /* check if all atoms in dihedral are either virtual sites, or used in
1205 construction of virtual sites. If so, keep it, if not throw away: */
1206 for(k=0; (k<4) && !bKeep; k++) { /* for all atoms in the dihedral */
1207 atom = ps->param[i].a[k];
1208 if (vsite_type[atom]==NOTSET) {
1210 for(m=0; (m<vsnral) && !bUsed; m++)
1211 if (atom == vsiteatoms[m])
1215 if (debug) fprintf(debug,"unused atom in dih: %u\n",atom+1);
1221 memcpy(&(ps->param[kept_i]),
1222 &(ps->param[i]),(size_t)sizeof(ps->param[0]));
1227 if (ps->nr != kept_i)
1228 fprintf(stderr,"Removed %4d %15ss with virtual sites, %5d left\n",
1229 ps->nr-kept_i, interaction_function[cftype].longname, kept_i);
1233 void clean_vsite_bondeds(t_params *plist, int natoms, gmx_bool bRmVSiteBds)
1235 int i,k,nvsite,ftype,vsite,parnr;
1238 at2vsitecon_t *at2vc;
1240 pindex=0; /* avoid warnings */
1241 /* make vsite_type array */
1242 snew(vsite_type,natoms);
1243 for(i=0; i<natoms; i++)
1244 vsite_type[i]=NOTSET;
1246 for(ftype=0; ftype<F_NRE; ftype++)
1247 if (interaction_function[ftype].flags & IF_VSITE) {
1248 nvsite+=plist[ftype].nr;
1250 while (i < plist[ftype].nr) {
1251 vsite = plist[ftype].param[i].AI;
1252 if ( vsite_type[vsite] == NOTSET)
1253 vsite_type[vsite] = ftype;
1255 gmx_fatal(FARGS,"multiple vsite constructions for atom %d",vsite+1);
1256 if (ftype == F_VSITEN) {
1257 while (i < plist[ftype].nr && plist[ftype].param[i].AI == vsite)
1265 /* the rest only if we have virtual sites: */
1267 fprintf(stderr,"Cleaning up constraints %swith virtual sites\n",
1268 bRmVSiteBds?"and constant bonded interactions ":"");
1270 /* Make a reverse list to avoid ninteractions^2 operations */
1271 at2vc = make_at2vsitecon(natoms,plist);
1273 snew(pindex,natoms);
1274 for(ftype=0; ftype<F_NRE; ftype++) {
1275 if ((interaction_function[ftype].flags & IF_VSITE) &&
1276 ftype != F_VSITEN) {
1277 for (parnr=0; (parnr<plist[ftype].nr); parnr++) {
1278 k=plist[ftype].param[parnr].AI;
1279 pindex[k].ftype=ftype;
1280 pindex[k].parnr=parnr;
1286 for(i=0; i<natoms; i++)
1287 fprintf(debug,"atom %d vsite_type %s\n",i,
1288 vsite_type[i]==NOTSET ? "NOTSET" :
1289 interaction_function[vsite_type[i]].name);
1291 /* remove things with vsite atoms */
1292 for(ftype=0; ftype<F_NRE; ftype++)
1293 if ( ( ( interaction_function[ftype].flags & IF_BOND ) && bRmVSiteBds ) ||
1294 ( interaction_function[ftype].flags & IF_CONSTRAINT ) ) {
1295 if (interaction_function[ftype].flags & (IF_BTYPE | IF_CONSTRAINT) )
1296 clean_vsite_bonds (plist, pindex, ftype, vsite_type);
1297 else if (interaction_function[ftype].flags & IF_ATYPE)
1298 clean_vsite_angles(plist, pindex, ftype, vsite_type, at2vc);
1299 else if ( (ftype==F_PDIHS) || (ftype==F_IDIHS) )
1300 clean_vsite_dihs (plist, pindex, ftype, vsite_type);
1302 /* check if we have constraints left with virtual sites in them */
1303 for(ftype=0; ftype<F_NRE; ftype++)
1304 if (interaction_function[ftype].flags & IF_CONSTRAINT)
1305 check_vsite_constraints(plist, ftype, vsite_type);
1307 done_at2vsitecon(natoms,at2vc);