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,2017, 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.
39 #include "vsite_parm.h"
48 #include "gromacs/gmxpreprocess/add_par.h"
49 #include "gromacs/gmxpreprocess/notset.h"
50 #include "gromacs/gmxpreprocess/resall.h"
51 #include "gromacs/gmxpreprocess/toputil.h"
52 #include "gromacs/math/functions.h"
53 #include "gromacs/math/units.h"
54 #include "gromacs/math/vec.h"
55 #include "gromacs/mdtypes/md_enums.h"
56 #include "gromacs/topology/ifunc.h"
57 #include "gromacs/topology/topology.h"
58 #include "gromacs/utility/cstringutil.h"
59 #include "gromacs/utility/fatalerror.h"
60 #include "gromacs/utility/gmxassert.h"
61 #include "gromacs/utility/smalloc.h"
62 #include "gromacs/utility/strconvert.h"
67 t_iatom
&ai() { return a
[0]; }
68 t_iatom
&aj() { return a
[1]; }
69 t_iatom
&ak() { return a
[2]; }
70 t_iatom
&al() { return a
[3]; }
81 vsitebondparam_t
*vsbp
;
89 static int vsite_bond_nrcheck(int ftype
)
93 if ((interaction_function
[ftype
].flags
& (IF_BTYPE
| IF_CONSTRAINT
| IF_ATYPE
)) || (ftype
== F_IDIHS
))
95 nrcheck
= NRAL(ftype
);
105 static void enter_bonded(int nratoms
, int *nrbonded
, t_mybonded
**bondeds
,
110 srenew(*bondeds
, *nrbonded
+1);
112 /* copy atom numbers */
113 for (j
= 0; j
< nratoms
; j
++)
115 (*bondeds
)[*nrbonded
].a
[j
] = param
->a
[j
];
118 (*bondeds
)[*nrbonded
].c
= param
->c0();
123 static void get_bondeds(int nrat
, t_iatom atoms
[],
124 at2vsitebond_t
*at2vb
,
125 int *nrbond
, t_mybonded
**bonds
,
126 int *nrang
, t_mybonded
**angles
,
127 int *nridih
, t_mybonded
**idihs
)
129 int k
, i
, ftype
, nrcheck
;
132 for (k
= 0; k
< nrat
; k
++)
134 for (i
= 0; i
< at2vb
[atoms
[k
]].nr
; i
++)
136 ftype
= at2vb
[atoms
[k
]].vsbp
[i
].ftype
;
137 param
= at2vb
[atoms
[k
]].vsbp
[i
].param
;
138 nrcheck
= vsite_bond_nrcheck(ftype
);
139 /* abuse nrcheck to see if we're adding bond, angle or idih */
142 case 2: enter_bonded(nrcheck
, nrbond
, bonds
, param
); break;
143 case 3: enter_bonded(nrcheck
, nrang
, angles
, param
); break;
144 case 4: enter_bonded(nrcheck
, nridih
, idihs
, param
); break;
150 static at2vsitebond_t
*make_at2vsitebond(int natoms
, t_params plist
[])
153 int ftype
, i
, j
, nrcheck
, nr
;
155 at2vsitebond_t
*at2vb
;
160 for (ftype
= 0; (ftype
< F_NRE
); ftype
++)
162 if ((interaction_function
[ftype
].flags
& IF_VSITE
) && ftype
!= F_VSITEN
)
164 for (i
= 0; (i
< plist
[ftype
].nr
); i
++)
166 for (j
= 0; j
< NRAL(ftype
); j
++)
168 bVSI
[plist
[ftype
].param
[i
].a
[j
]] = TRUE
;
174 for (ftype
= 0; (ftype
< F_NRE
); ftype
++)
176 nrcheck
= vsite_bond_nrcheck(ftype
);
179 for (i
= 0; (i
< plist
[ftype
].nr
); i
++)
181 aa
= plist
[ftype
].param
[i
].a
;
182 for (j
= 0; j
< nrcheck
; j
++)
186 nr
= at2vb
[aa
[j
]].nr
;
189 srenew(at2vb
[aa
[j
]].vsbp
, nr
+10);
191 at2vb
[aa
[j
]].vsbp
[nr
].ftype
= ftype
;
192 at2vb
[aa
[j
]].vsbp
[nr
].param
= &plist
[ftype
].param
[i
];
204 static void done_at2vsitebond(int natoms
, at2vsitebond_t
*at2vb
)
208 for (i
= 0; i
< natoms
; i
++)
212 sfree(at2vb
[i
].vsbp
);
218 static at2vsitecon_t
*make_at2vsitecon(int natoms
, t_params plist
[])
221 int ftype
, i
, j
, ai
, aj
, nr
;
222 at2vsitecon_t
*at2vc
;
227 for (ftype
= 0; (ftype
< F_NRE
); ftype
++)
229 if ((interaction_function
[ftype
].flags
& IF_VSITE
) && ftype
!= F_VSITEN
)
231 for (i
= 0; (i
< plist
[ftype
].nr
); i
++)
233 for (j
= 0; j
< NRAL(ftype
); j
++)
235 bVSI
[plist
[ftype
].param
[i
].a
[j
]] = TRUE
;
241 for (ftype
= 0; (ftype
< F_NRE
); ftype
++)
243 if (interaction_function
[ftype
].flags
& IF_CONSTRAINT
)
245 for (i
= 0; (i
< plist
[ftype
].nr
); i
++)
247 ai
= plist
[ftype
].param
[i
].ai();
248 aj
= plist
[ftype
].param
[i
].aj();
249 if (bVSI
[ai
] && bVSI
[aj
])
251 /* Store forward direction */
255 srenew(at2vc
[ai
].aj
, nr
+10);
257 at2vc
[ai
].aj
[nr
] = aj
;
259 /* Store backward direction */
263 srenew(at2vc
[aj
].aj
, nr
+10);
265 at2vc
[aj
].aj
[nr
] = ai
;
276 static void done_at2vsitecon(int natoms
, at2vsitecon_t
*at2vc
)
280 for (i
= 0; i
< natoms
; i
++)
291 static void print_bad(FILE *fp
,
292 int nrbond
, t_mybonded
*bonds
,
293 int nrang
, t_mybonded
*angles
,
294 int nridih
, t_mybonded
*idihs
)
300 fprintf(fp
, "bonds:");
301 for (i
= 0; i
< nrbond
; i
++)
303 fprintf(fp
, " %d-%d (%g)",
304 bonds
[i
].ai()+1, bonds
[i
].aj()+1, bonds
[i
].c
);
310 fprintf(fp
, "angles:");
311 for (i
= 0; i
< nrang
; i
++)
313 fprintf(fp
, " %d-%d-%d (%g)",
314 angles
[i
].ai()+1, angles
[i
].aj()+1,
315 angles
[i
].ak()+1, angles
[i
].c
);
321 fprintf(fp
, "idihs:");
322 for (i
= 0; i
< nridih
; i
++)
324 fprintf(fp
, " %d-%d-%d-%d (%g)",
325 idihs
[i
].ai()+1, idihs
[i
].aj()+1,
326 idihs
[i
].ak()+1, idihs
[i
].al()+1, idihs
[i
].c
);
332 static void print_param(FILE *fp
, int ftype
, int i
, t_param
*param
)
335 static int prev_ftype
= NOTSET
;
336 static int prev_i
= NOTSET
;
339 if ( (ftype
!= prev_ftype
) || (i
!= prev_i
) )
345 fprintf(fp
, "(%d) plist[%s].param[%d]",
346 pass
, interaction_function
[ftype
].name
, i
);
347 for (j
= 0; j
< NRFP(ftype
); j
++)
349 fprintf(fp
, ".c[%d]=%g ", j
, param
->c
[j
]);
355 static real
get_bond_length(int nrbond
, t_mybonded bonds
[],
356 t_iatom ai
, t_iatom aj
)
362 for (i
= 0; i
< nrbond
&& (bondlen
== NOTSET
); i
++)
364 /* check both ways */
365 if ( ( (ai
== bonds
[i
].ai()) && (aj
== bonds
[i
].aj()) ) ||
366 ( (ai
== bonds
[i
].aj()) && (aj
== bonds
[i
].ai()) ) )
368 bondlen
= bonds
[i
].c
; /* note: bonds[i].c might be NOTSET */
374 static real
get_angle(int nrang
, t_mybonded angles
[],
375 t_iatom ai
, t_iatom aj
, t_iatom ak
)
381 for (i
= 0; i
< nrang
&& (angle
== NOTSET
); i
++)
383 /* check both ways */
384 if ( ( (ai
== angles
[i
].ai()) && (aj
== angles
[i
].aj()) && (ak
== angles
[i
].ak()) ) ||
385 ( (ai
== angles
[i
].ak()) && (aj
== angles
[i
].aj()) && (ak
== angles
[i
].ai()) ) )
387 angle
= DEG2RAD
*angles
[i
].c
;
393 static char *get_atomtype_name_AB(t_atom
*atom
, gpp_atomtype_t atype
)
397 name
= get_atomtype_name(atom
->type
, atype
);
399 /* When using the decoupling option, atom types are changed
400 * to decoupled for the non-bonded interactions, but the virtual
401 * sites constructions should be based on the "normal" interactions.
402 * So we return the state B atom type name if the state A atom
403 * type is the decoupled one. We should actually check for the atom
404 * type number, but that's not passed here. So we check for
405 * the decoupled atom type name. This should not cause trouble
406 * as this code is only used for topologies with v-sites without
407 * parameters generated by pdb2gmx.
409 if (strcmp(name
, "decoupled") == 0)
411 name
= get_atomtype_name(atom
->typeB
, atype
);
417 static gmx_bool
calc_vsite3_param(gpp_atomtype_t atype
,
418 t_param
*param
, t_atoms
*at
,
419 int nrbond
, t_mybonded
*bonds
,
420 int nrang
, t_mybonded
*angles
)
422 /* i = virtual site | ,k
423 * j = 1st bonded heavy atom | i-j
424 * k,l = 2nd bonded atoms | `l
427 gmx_bool bXH3
, bError
;
428 real bjk
, bjl
, a
= -1, b
= -1;
429 /* check if this is part of a NH3 , NH2-umbrella or CH3 group,
430 * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
434 for (i
= 0; i
< 4; i
++)
436 fprintf(debug
, "atom %d type %s ",
438 get_atomtype_name_AB(&at
->atom
[param
->a
[i
]], atype
));
440 fprintf(debug
, "\n");
443 ( (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->ak()], atype
), "MNH", 3) == 0) &&
444 (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->al()], atype
), "MNH", 3) == 0) ) ||
445 ( (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->ak()], atype
), "MCH3", 4) == 0) &&
446 (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->al()], atype
), "MCH3", 4) == 0) );
448 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
449 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
450 bError
= (bjk
== NOTSET
) || (bjl
== NOTSET
);
453 /* now we get some XH2/XH3 group specific construction */
454 /* note: we call the heavy atom 'C' and the X atom 'N' */
455 real bMM
, bCM
, bCN
, bNH
, aCNH
, dH
, rH
, dM
, rM
;
458 /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
459 bError
= bError
|| (bjk
!= bjl
);
461 /* the X atom (C or N) in the XH2/XH3 group is the first after the masses: */
462 aN
= std::max(param
->ak(), param
->al())+1;
464 /* get common bonds */
465 bMM
= get_bond_length(nrbond
, bonds
, param
->ak(), param
->al());
467 bCN
= get_bond_length(nrbond
, bonds
, param
->aj(), aN
);
468 bError
= bError
|| (bMM
== NOTSET
) || (bCN
== NOTSET
);
470 /* calculate common things */
472 dM
= std::sqrt( gmx::square(bCM
) - gmx::square(rM
) );
474 /* are we dealing with the X atom? */
475 if (param
->ai() == aN
)
477 /* this is trivial */
478 a
= b
= 0.5 * bCN
/dM
;
483 /* get other bondlengths and angles: */
484 bNH
= get_bond_length(nrbond
, bonds
, aN
, param
->ai());
485 aCNH
= get_angle (nrang
, angles
, param
->aj(), aN
, param
->ai());
486 bError
= bError
|| (bNH
== NOTSET
) || (aCNH
== NOTSET
);
489 dH
= bCN
- bNH
* std::cos(aCNH
);
490 rH
= bNH
* std::sin(aCNH
);
492 a
= 0.5 * ( dH
/dM
+ rH
/rM
);
493 b
= 0.5 * ( dH
/dM
- rH
/rM
);
498 gmx_fatal(FARGS
, "calc_vsite3_param not implemented for the general case "
499 "(atom %d)", param
->ai()+1);
507 fprintf(debug
, "params for vsite3 %d: %g %g\n",
508 param
->ai()+1, param
->c0(), param
->c1());
514 static gmx_bool
calc_vsite3fd_param(t_param
*param
,
515 int nrbond
, t_mybonded
*bonds
,
516 int nrang
, t_mybonded
*angles
)
518 /* i = virtual site | ,k
519 * j = 1st bonded heavy atom | i-j
520 * k,l = 2nd bonded atoms | `l
524 real bij
, bjk
, bjl
, aijk
, aijl
, rk
, rl
;
526 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
527 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
528 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
529 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
530 aijl
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->al());
531 bError
= (bij
== NOTSET
) || (bjk
== NOTSET
) || (bjl
== NOTSET
) ||
532 (aijk
== NOTSET
) || (aijl
== NOTSET
);
534 rk
= bjk
* std::sin(aijk
);
535 rl
= bjl
* std::sin(aijl
);
536 param
->c0() = rk
/ (rk
+ rl
);
537 param
->c1() = -bij
; /* 'bond'-length for fixed distance vsite */
541 fprintf(debug
, "params for vsite3fd %d: %g %g\n",
542 param
->ai()+1, param
->c0(), param
->c1());
547 static gmx_bool
calc_vsite3fad_param(t_param
*param
,
548 int nrbond
, t_mybonded
*bonds
,
549 int nrang
, t_mybonded
*angles
)
551 /* i = virtual site |
552 * j = 1st bonded heavy atom | i-j
553 * k = 2nd bonded heavy atom | `k-l
554 * l = 3d bonded heavy atom |
557 gmx_bool bSwapParity
, bError
;
560 bSwapParity
= ( param
->c1() == -1 );
562 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
563 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
564 bError
= (bij
== NOTSET
) || (aijk
== NOTSET
);
566 param
->c1() = bij
; /* 'bond'-length for fixed distance vsite */
567 param
->c0() = RAD2DEG
*aijk
; /* 'bond'-angle for fixed angle vsite */
571 param
->c0() = 360 - param
->c0();
576 fprintf(debug
, "params for vsite3fad %d: %g %g\n",
577 param
->ai()+1, param
->c0(), param
->c1());
582 static gmx_bool
calc_vsite3out_param(gpp_atomtype_t atype
,
583 t_param
*param
, t_atoms
*at
,
584 int nrbond
, t_mybonded
*bonds
,
585 int nrang
, t_mybonded
*angles
)
587 /* i = virtual site | ,k
588 * j = 1st bonded heavy atom | i-j
589 * k,l = 2nd bonded atoms | `l
590 * NOTE: i is out of the j-k-l plane!
593 gmx_bool bXH3
, bError
, bSwapParity
;
594 real bij
, bjk
, bjl
, aijk
, aijl
, akjl
, pijk
, pijl
, a
, b
, c
;
596 /* check if this is part of a NH2-umbrella, NH3 or CH3 group,
597 * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
601 for (i
= 0; i
< 4; i
++)
603 fprintf(debug
, "atom %d type %s ",
604 param
->a
[i
]+1, get_atomtype_name_AB(&at
->atom
[param
->a
[i
]], atype
));
606 fprintf(debug
, "\n");
609 ( (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->ak()], atype
), "MNH", 3) == 0) &&
610 (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->al()], atype
), "MNH", 3) == 0) ) ||
611 ( (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->ak()], atype
), "MCH3", 4) == 0) &&
612 (gmx_strncasecmp(get_atomtype_name_AB(&at
->atom
[param
->al()], atype
), "MCH3", 4) == 0) );
614 /* check if construction parity must be swapped */
615 bSwapParity
= ( param
->c1() == -1 );
617 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
618 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
619 bError
= (bjk
== NOTSET
) || (bjl
== NOTSET
);
622 /* now we get some XH3 group specific construction */
623 /* note: we call the heavy atom 'C' and the X atom 'N' */
624 real bMM
, bCM
, bCN
, bNH
, aCNH
, dH
, rH
, rHx
, rHy
, dM
, rM
;
627 /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
628 bError
= bError
|| (bjk
!= bjl
);
630 /* the X atom (C or N) in the XH3 group is the first after the masses: */
631 aN
= std::max(param
->ak(), param
->al())+1;
633 /* get all bondlengths and angles: */
634 bMM
= get_bond_length(nrbond
, bonds
, param
->ak(), param
->al());
636 bCN
= get_bond_length(nrbond
, bonds
, param
->aj(), aN
);
637 bNH
= get_bond_length(nrbond
, bonds
, aN
, param
->ai());
638 aCNH
= get_angle (nrang
, angles
, param
->aj(), aN
, param
->ai());
640 (bMM
== NOTSET
) || (bCN
== NOTSET
) || (bNH
== NOTSET
) || (aCNH
== NOTSET
);
643 dH
= bCN
- bNH
* std::cos(aCNH
);
644 rH
= bNH
* std::sin(aCNH
);
645 /* we assume the H's are symmetrically distributed */
646 rHx
= rH
*std::cos(DEG2RAD
*30);
647 rHy
= rH
*std::sin(DEG2RAD
*30);
649 dM
= std::sqrt( gmx::square(bCM
) - gmx::square(rM
) );
650 a
= 0.5*( (dH
/dM
) - (rHy
/rM
) );
651 b
= 0.5*( (dH
/dM
) + (rHy
/rM
) );
657 /* this is the general construction */
659 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
660 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
661 aijl
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->al());
662 akjl
= get_angle (nrang
, angles
, param
->ak(), param
->aj(), param
->al());
664 (bij
== NOTSET
) || (aijk
== NOTSET
) || (aijl
== NOTSET
) || (akjl
== NOTSET
);
666 pijk
= std::cos(aijk
)*bij
;
667 pijl
= std::cos(aijl
)*bij
;
668 a
= ( pijk
+ (pijk
*std::cos(akjl
)-pijl
) * std::cos(akjl
) / gmx::square(std::sin(akjl
)) ) / bjk
;
669 b
= ( pijl
+ (pijl
*std::cos(akjl
)-pijk
) * std::cos(akjl
) / gmx::square(std::sin(akjl
)) ) / bjl
;
670 c
= -std::sqrt( gmx::square(bij
) -
671 ( gmx::square(pijk
) - 2*pijk
*pijl
*std::cos(akjl
) + gmx::square(pijl
) )
672 / gmx::square(std::sin(akjl
)) )
673 / ( bjk
*bjl
*std::sin(akjl
) );
688 fprintf(debug
, "params for vsite3out %d: %g %g %g\n",
689 param
->ai()+1, param
->c0(), param
->c1(), param
->c2());
694 static gmx_bool
calc_vsite4fd_param(t_param
*param
,
695 int nrbond
, t_mybonded
*bonds
,
696 int nrang
, t_mybonded
*angles
)
698 /* i = virtual site | ,k
699 * j = 1st bonded heavy atom | i-j-m
700 * k,l,m = 2nd bonded atoms | `l
704 real bij
, bjk
, bjl
, bjm
, aijk
, aijl
, aijm
, akjm
, akjl
;
705 real pk
, pl
, pm
, cosakl
, cosakm
, sinakl
, sinakm
, cl
, cm
;
707 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
708 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
709 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
710 bjm
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->am());
711 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
712 aijl
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->al());
713 aijm
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->am());
714 akjm
= get_angle (nrang
, angles
, param
->ak(), param
->aj(), param
->am());
715 akjl
= get_angle (nrang
, angles
, param
->ak(), param
->aj(), param
->al());
716 bError
= (bij
== NOTSET
) || (bjk
== NOTSET
) || (bjl
== NOTSET
) || (bjm
== NOTSET
) ||
717 (aijk
== NOTSET
) || (aijl
== NOTSET
) || (aijm
== NOTSET
) || (akjm
== NOTSET
) ||
722 pk
= bjk
*std::sin(aijk
);
723 pl
= bjl
*std::sin(aijl
);
724 pm
= bjm
*std::sin(aijm
);
725 cosakl
= (std::cos(akjl
) - std::cos(aijk
)*std::cos(aijl
)) / (std::sin(aijk
)*std::sin(aijl
));
726 cosakm
= (std::cos(akjm
) - std::cos(aijk
)*std::cos(aijm
)) / (std::sin(aijk
)*std::sin(aijm
));
727 if (cosakl
< -1 || cosakl
> 1 || cosakm
< -1 || cosakm
> 1)
729 fprintf(stderr
, "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
730 param
->ai()+1, RAD2DEG
*aijk
, RAD2DEG
*aijl
, RAD2DEG
*aijm
);
731 gmx_fatal(FARGS
, "invalid construction in calc_vsite4fd for atom %d: "
732 "cosakl=%f, cosakm=%f\n", param
->ai()+1, cosakl
, cosakm
);
734 sinakl
= std::sqrt(1-gmx::square(cosakl
));
735 sinakm
= std::sqrt(1-gmx::square(cosakm
));
737 /* note: there is a '+' because of the way the sines are calculated */
738 cl
= -pk
/ ( pl
*cosakl
- pk
+ pl
*sinakl
*(pm
*cosakm
-pk
)/(pm
*sinakm
) );
739 cm
= -pk
/ ( pm
*cosakm
- pk
+ pm
*sinakm
*(pl
*cosakl
-pk
)/(pl
*sinakl
) );
746 fprintf(debug
, "params for vsite4fd %d: %g %g %g\n",
747 param
->ai()+1, param
->c0(), param
->c1(), param
->c2());
756 calc_vsite4fdn_param(t_param
*param
,
757 int nrbond
, t_mybonded
*bonds
,
758 int nrang
, t_mybonded
*angles
)
760 /* i = virtual site | ,k
761 * j = 1st bonded heavy atom | i-j-m
762 * k,l,m = 2nd bonded atoms | `l
766 real bij
, bjk
, bjl
, bjm
, aijk
, aijl
, aijm
;
767 real pk
, pl
, pm
, a
, b
;
769 bij
= get_bond_length(nrbond
, bonds
, param
->ai(), param
->aj());
770 bjk
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->ak());
771 bjl
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->al());
772 bjm
= get_bond_length(nrbond
, bonds
, param
->aj(), param
->am());
773 aijk
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->ak());
774 aijl
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->al());
775 aijm
= get_angle (nrang
, angles
, param
->ai(), param
->aj(), param
->am());
777 bError
= (bij
== NOTSET
) || (bjk
== NOTSET
) || (bjl
== NOTSET
) || (bjm
== NOTSET
) ||
778 (aijk
== NOTSET
) || (aijl
== NOTSET
) || (aijm
== NOTSET
);
783 /* Calculate component of bond j-k along the direction i-j */
784 pk
= -bjk
*std::cos(aijk
);
786 /* Calculate component of bond j-l along the direction i-j */
787 pl
= -bjl
*std::cos(aijl
);
789 /* Calculate component of bond j-m along the direction i-j */
790 pm
= -bjm
*std::cos(aijm
);
792 if (fabs(pl
) < 1000*GMX_REAL_MIN
|| fabs(pm
) < 1000*GMX_REAL_MIN
)
794 fprintf(stderr
, "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f\n",
795 param
->ai()+1, RAD2DEG
*aijk
, RAD2DEG
*aijl
, RAD2DEG
*aijm
);
796 gmx_fatal(FARGS
, "invalid construction in calc_vsite4fdn for atom %d: "
797 "pl=%f, pm=%f\n", param
->ai()+1, pl
, pm
);
809 fprintf(debug
, "params for vsite4fdn %d: %g %g %g\n",
810 param
->ai()+1, param
->c0(), param
->c1(), param
->c2());
819 int set_vsites(gmx_bool bVerbose
, t_atoms
*atoms
, gpp_atomtype_t atype
,
823 int nvsite
, nrbond
, nrang
, nridih
, nrset
;
824 gmx_bool bFirst
, bSet
, bERROR
;
825 at2vsitebond_t
*at2vb
;
834 fprintf(debug
, "\nCalculating parameters for virtual sites\n");
837 /* Make a reverse list to avoid ninteractions^2 operations */
838 at2vb
= make_at2vsitebond(atoms
->nr
, plist
);
840 for (ftype
= 0; (ftype
< F_NRE
); ftype
++)
842 if (interaction_function
[ftype
].flags
& IF_VSITE
)
844 nvsite
+= plist
[ftype
].nr
;
846 if (ftype
== F_VSITEN
)
848 /* We don't calculate parameters for VSITEN */
853 for (i
= 0; (i
< plist
[ftype
].nr
); i
++)
855 /* check if all parameters are set */
857 for (j
= 0; j
< NRFP(ftype
) && bSet
; j
++)
859 bSet
= plist
[ftype
].param
[i
].c
[j
] != NOTSET
;
864 fprintf(debug
, "bSet=%s ", gmx::boolToString(bSet
));
865 print_param(debug
, ftype
, i
, &plist
[ftype
].param
[i
]);
869 if (bVerbose
&& bFirst
)
871 fprintf(stderr
, "Calculating parameters for virtual sites\n");
875 nrbond
= nrang
= nridih
= 0;
880 /* now set the vsite parameters: */
881 get_bondeds(NRAL(ftype
), plist
[ftype
].param
[i
].a
, at2vb
,
882 &nrbond
, &bonds
, &nrang
, &angles
, &nridih
, &idihs
);
885 fprintf(debug
, "Found %d bonds, %d angles and %d idihs "
886 "for virtual site %d (%s)\n", nrbond
, nrang
, nridih
,
887 plist
[ftype
].param
[i
].ai()+1,
888 interaction_function
[ftype
].longname
);
889 print_bad(debug
, nrbond
, bonds
, nrang
, angles
, nridih
, idihs
);
895 calc_vsite3_param(atype
, &(plist
[ftype
].param
[i
]), atoms
,
896 nrbond
, bonds
, nrang
, angles
);
900 calc_vsite3fd_param(&(plist
[ftype
].param
[i
]),
901 nrbond
, bonds
, nrang
, angles
);
905 calc_vsite3fad_param(&(plist
[ftype
].param
[i
]),
906 nrbond
, bonds
, nrang
, angles
);
910 calc_vsite3out_param(atype
, &(plist
[ftype
].param
[i
]), atoms
,
911 nrbond
, bonds
, nrang
, angles
);
915 calc_vsite4fd_param(&(plist
[ftype
].param
[i
]),
916 nrbond
, bonds
, nrang
, angles
);
920 calc_vsite4fdn_param(&(plist
[ftype
].param
[i
]),
921 nrbond
, bonds
, nrang
, angles
);
924 gmx_fatal(FARGS
, "Automatic parameter generation not supported "
926 interaction_function
[ftype
].longname
,
927 plist
[ftype
].param
[i
].ai()+1);
932 gmx_fatal(FARGS
, "Automatic parameter generation not supported "
933 "for %s atom %d for this bonding configuration",
934 interaction_function
[ftype
].longname
,
935 plist
[ftype
].param
[i
].ai()+1);
942 if (debug
&& plist
[ftype
].nr
)
944 fprintf(stderr
, "Calculated parameters for %d out of %d %s atoms\n",
945 nrset
, plist
[ftype
].nr
, interaction_function
[ftype
].longname
);
950 done_at2vsitebond(atoms
->nr
, at2vb
);
955 void set_vsites_ptype(gmx_bool bVerbose
, gmx_moltype_t
*molt
)
964 fprintf(stderr
, "Setting particle type to V for virtual sites\n");
968 fprintf(stderr
, "checking %d functypes\n", F_NRE
);
970 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
972 il
= &molt
->ilist
[ftype
];
973 if (interaction_function
[ftype
].flags
& IF_VSITE
)
975 nra
= interaction_function
[ftype
].nratoms
;
981 fprintf(stderr
, "doing %d %s virtual sites\n",
982 (nrd
/ (nra
+1)), interaction_function
[ftype
].longname
);
985 for (i
= 0; (i
< nrd
); )
987 /* The virtual site */
989 molt
->atoms
.atom
[avsite
].ptype
= eptVSite
;
1003 static void check_vsite_constraints(t_params
*plist
,
1004 int cftype
, int vsite_type
[])
1011 ps
= &(plist
[cftype
]);
1012 for (i
= 0; (i
< ps
->nr
); i
++)
1014 for (k
= 0; k
< 2; k
++)
1016 atom
= ps
->param
[i
].a
[k
];
1017 if (vsite_type
[atom
] != NOTSET
)
1019 fprintf(stderr
, "ERROR: Cannot have constraint (%d-%d) with virtual site (%d)\n",
1020 ps
->param
[i
].ai()+1, ps
->param
[i
].aj()+1, atom
+1);
1027 gmx_fatal(FARGS
, "There were %d virtual sites involved in constraints", n
);
1031 static void clean_vsite_bonds(t_params
*plist
, t_pindex pindex
[],
1032 int cftype
, int vsite_type
[])
1034 int ftype
, i
, j
, k
, m
, n
, nvsite
, nOut
, kept_i
;
1035 int nconverted
, nremoved
;
1036 int atom
, oatom
, at1
, at2
;
1037 gmx_bool bKeep
, bRemove
, bUsed
, bPresent
, bThisFD
, bThisOUT
, bAllFD
, bFirstTwo
;
1040 if (cftype
== F_CONNBONDS
)
1045 ps
= &(plist
[cftype
]);
1050 for (i
= 0; (i
< ps
->nr
); i
++) /* for all bonds in the plist */
1053 const int *first_atoms
= nullptr;
1058 /* check if all virtual sites are constructed from the same atoms */
1062 fprintf(debug
, "constr %d %d:", ps
->param
[i
].ai()+1, ps
->param
[i
].aj()+1);
1064 for (k
= 0; (k
< 2) && !bKeep
&& !bRemove
; k
++)
1066 /* for all atoms in the bond */
1067 atom
= ps
->param
[i
].a
[k
];
1068 if (vsite_type
[atom
] != NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1071 bThisFD
= ( (pindex
[atom
].ftype
== F_VSITE3FD
) ||
1072 (pindex
[atom
].ftype
== F_VSITE3FAD
) ||
1073 (pindex
[atom
].ftype
== F_VSITE4FD
) ||
1074 (pindex
[atom
].ftype
== F_VSITE4FDN
) );
1075 bThisOUT
= ( (pindex
[atom
].ftype
== F_VSITE3OUT
) &&
1076 (interaction_function
[cftype
].flags
& IF_CONSTRAINT
) );
1077 bAllFD
= bAllFD
&& bThisFD
;
1078 if (bThisFD
|| bThisOUT
)
1082 fprintf(debug
, " %s", bThisOUT
? "out" : "fd");
1084 oatom
= ps
->param
[i
].a
[1-k
]; /* the other atom */
1085 if (vsite_type
[oatom
] == NOTSET
&&
1086 vsite_type
[oatom
] != F_VSITEN
&&
1087 oatom
== plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].aj())
1089 /* if the other atom isn't a vsite, and it is AI */
1097 fprintf(debug
, " D-AI");
1103 /* TODO This fragment, and corresponding logic in
1104 clean_vsite_angles and clean_vsite_dihs should
1105 be refactored into a common function */
1108 /* if this is the first vsite we encounter then
1109 store construction atoms */
1110 /* TODO This would be nicer to implement with
1111 a C++ "vector view" class" with an
1112 STL-container-like interface. */
1113 vsnral
= NRAL(pindex
[atom
].ftype
) - 1;
1114 first_atoms
= plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].a
+ 1;
1118 GMX_ASSERT(vsnral
!= 0, "nvsite > 1 must have vsnral != 0");
1119 GMX_ASSERT(first_atoms
!= NULL
, "nvsite > 1 must have first_atoms != NULL");
1120 /* if it is not the first then
1121 check if this vsite is constructed from the same atoms */
1122 if (vsnral
== NRAL(pindex
[atom
].ftype
)-1)
1124 for (m
= 0; (m
< vsnral
) && !bKeep
; m
++)
1129 atoms
= plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].a
+ 1;
1130 for (n
= 0; (n
< vsnral
) && !bPresent
; n
++)
1132 if (atoms
[m
] == first_atoms
[n
])
1142 fprintf(debug
, " !present");
1152 fprintf(debug
, " !same#at");
1166 /* if we have no virtual sites in this bond, keep it */
1171 fprintf(debug
, " no vsite");
1176 /* TODO This loop and the corresponding loop in
1177 check_vsite_angles should be refactored into a common
1179 /* check if all non-vsite atoms are used in construction: */
1181 for (k
= 0; (k
< 2) && !bKeep
; k
++) /* for all atoms in the bond */
1183 atom
= ps
->param
[i
].a
[k
];
1184 if (vsite_type
[atom
] == NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1187 for (m
= 0; (m
< vsnral
) && !bUsed
; m
++)
1189 GMX_ASSERT(first_atoms
!= NULL
, "If we've seen a vsite before, we know what its first atom index was");
1191 if (atom
== first_atoms
[m
])
1194 bFirstTwo
= bFirstTwo
&& m
< 2;
1202 fprintf(debug
, " !used");
1208 if (!( bAllFD
&& bFirstTwo
) )
1210 /* Two atom bonded interactions include constraints.
1211 * We need to remove constraints between vsite pairs that have
1212 * a fixed distance due to being constructed from the same
1213 * atoms, since this can be numerically unstable.
1215 for (m
= 0; m
< vsnral
&& !bKeep
; m
++) /* all constr. atoms */
1217 at1
= first_atoms
[m
];
1218 at2
= first_atoms
[(m
+1) % vsnral
];
1220 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1222 if (interaction_function
[ftype
].flags
& IF_CONSTRAINT
)
1224 for (j
= 0; (j
< plist
[ftype
].nr
) && !bPresent
; j
++)
1226 /* all constraints until one matches */
1227 bPresent
= ( ( (plist
[ftype
].param
[j
].ai() == at1
) &&
1228 (plist
[ftype
].param
[j
].aj() == at2
) ) ||
1229 ( (plist
[ftype
].param
[j
].ai() == at2
) &&
1230 (plist
[ftype
].param
[j
].aj() == at1
) ) );
1239 fprintf(debug
, " !bonded");
1250 fprintf(debug
, " keeping");
1252 /* now copy the bond to the new array */
1253 ps
->param
[kept_i
] = ps
->param
[i
];
1256 else if (IS_CHEMBOND(cftype
))
1258 srenew(plist
[F_CONNBONDS
].param
, plist
[F_CONNBONDS
].nr
+1);
1259 plist
[F_CONNBONDS
].param
[plist
[F_CONNBONDS
].nr
] = ps
->param
[i
];
1260 plist
[F_CONNBONDS
].nr
++;
1269 fprintf(debug
, "\n");
1275 fprintf(stderr
, "Removed %4d %15ss with virtual sites, %5d left\n",
1276 nremoved
, interaction_function
[cftype
].longname
, kept_i
);
1280 fprintf(stderr
, "Converted %4d %15ss with virtual sites to connections, %5d left\n",
1281 nconverted
, interaction_function
[cftype
].longname
, kept_i
);
1285 fprintf(stderr
, "Warning: removed %d %ss with vsite with %s construction\n"
1286 " This vsite construction does not guarantee constant "
1288 " If the constructions were generated by pdb2gmx ignore "
1290 nOut
, interaction_function
[cftype
].longname
,
1291 interaction_function
[F_VSITE3OUT
].longname
);
1296 static void clean_vsite_angles(t_params
*plist
, t_pindex pindex
[],
1297 int cftype
, int vsite_type
[],
1298 at2vsitecon_t
*at2vc
)
1300 int i
, j
, k
, m
, n
, nvsite
, kept_i
;
1302 gmx_bool bKeep
, bUsed
, bPresent
, bAll3FAD
, bFirstTwo
;
1305 ps
= &(plist
[cftype
]);
1307 for (i
= 0; (i
< ps
->nr
); i
++) /* for all angles in the plist */
1310 const int *first_atoms
= nullptr;
1314 /* check if all virtual sites are constructed from the same atoms */
1316 for (k
= 0; (k
< 3) && !bKeep
; k
++) /* for all atoms in the angle */
1318 atom
= ps
->param
[i
].a
[k
];
1319 if (vsite_type
[atom
] != NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1322 bAll3FAD
= bAll3FAD
&& (pindex
[atom
].ftype
== F_VSITE3FAD
);
1325 /* store construction atoms of first vsite */
1326 vsnral
= NRAL(pindex
[atom
].ftype
) - 1;
1327 first_atoms
= plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].a
+ 1;
1331 GMX_ASSERT(vsnral
!= 0, "If we've seen a vsite before, we know how many constructing atoms it had");
1332 GMX_ASSERT(first_atoms
!= NULL
, "If we've seen a vsite before, we know what its first atom index was");
1333 /* check if this vsite is constructed from the same atoms */
1334 if (vsnral
== NRAL(pindex
[atom
].ftype
)-1)
1336 for (m
= 0; (m
< vsnral
) && !bKeep
; m
++)
1341 atoms
= plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].a
+ 1;
1342 for (n
= 0; (n
< vsnral
) && !bPresent
; n
++)
1344 if (atoms
[m
] == first_atoms
[n
])
1363 /* keep all angles with no virtual sites in them or
1364 with virtual sites with more than 3 constr. atoms */
1365 if (nvsite
== 0 && vsnral
> 3)
1370 /* check if all non-vsite atoms are used in construction: */
1372 for (k
= 0; (k
< 3) && !bKeep
; k
++) /* for all atoms in the angle */
1374 atom
= ps
->param
[i
].a
[k
];
1375 if (vsite_type
[atom
] == NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1378 for (m
= 0; (m
< vsnral
) && !bUsed
; m
++)
1380 GMX_ASSERT(first_atoms
!= NULL
, "If we've seen a vsite before, we know what its first atom index was");
1382 if (atom
== first_atoms
[m
])
1385 bFirstTwo
= bFirstTwo
&& m
< 2;
1395 if (!( bAll3FAD
&& bFirstTwo
) )
1397 /* check if all constructing atoms are constrained together */
1398 for (m
= 0; m
< vsnral
&& !bKeep
; m
++) /* all constr. atoms */
1400 at1
= first_atoms
[m
];
1401 at2
= first_atoms
[(m
+1) % vsnral
];
1403 for (j
= 0; j
< at2vc
[at1
].nr
; j
++)
1405 if (at2vc
[at1
].aj
[j
] == at2
)
1419 /* now copy the angle to the new array */
1420 ps
->param
[kept_i
] = ps
->param
[i
];
1425 if (ps
->nr
!= kept_i
)
1427 fprintf(stderr
, "Removed %4d %15ss with virtual sites, %5d left\n",
1428 ps
->nr
-kept_i
, interaction_function
[cftype
].longname
, kept_i
);
1433 static void clean_vsite_dihs(t_params
*plist
, t_pindex pindex
[],
1434 int cftype
, int vsite_type
[])
1439 ps
= &(plist
[cftype
]);
1442 for (i
= 0; (i
< ps
->nr
); i
++) /* for all dihedrals in the plist */
1444 int k
, m
, n
, nvsite
;
1446 const int *first_atoms
= nullptr;
1448 gmx_bool bKeep
, bUsed
, bPresent
;
1452 /* check if all virtual sites are constructed from the same atoms */
1454 for (k
= 0; (k
< 4) && !bKeep
; k
++) /* for all atoms in the dihedral */
1456 atom
= ps
->param
[i
].a
[k
];
1457 if (vsite_type
[atom
] != NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1461 /* store construction atoms of first vsite */
1462 vsnral
= NRAL(pindex
[atom
].ftype
) - 1;
1463 first_atoms
= plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].a
+ 1;
1466 fprintf(debug
, "dih w. vsite: %d %d %d %d\n",
1467 ps
->param
[i
].ai()+1, ps
->param
[i
].aj()+1,
1468 ps
->param
[i
].ak()+1, ps
->param
[i
].al()+1);
1469 fprintf(debug
, "vsite %d from: %d %d %d\n",
1470 atom
+1, first_atoms
[0]+1, first_atoms
[1]+1, first_atoms
[2]+1);
1475 GMX_ASSERT(vsnral
!= 0, "If we've seen a vsite before, we know how many constructing atoms it had");
1476 GMX_ASSERT(first_atoms
!= NULL
, "If we've seen a vsite before, we know what its first atom index was");
1477 /* check if this vsite is constructed from the same atoms */
1478 if (vsnral
== NRAL(pindex
[atom
].ftype
)-1)
1480 for (m
= 0; (m
< vsnral
) && !bKeep
; m
++)
1485 atoms
= plist
[pindex
[atom
].ftype
].param
[pindex
[atom
].parnr
].a
+ 1;
1486 for (n
= 0; (n
< vsnral
) && !bPresent
; n
++)
1488 if (atoms
[m
] == first_atoms
[n
])
1500 /* TODO clean_site_bonds and _angles do this increment
1501 at the top of the loop. Refactor this for
1507 /* keep all dihedrals with no virtual sites in them */
1513 /* check if all atoms in dihedral are either virtual sites, or used in
1514 construction of virtual sites. If so, keep it, if not throw away: */
1515 for (k
= 0; (k
< 4) && !bKeep
; k
++) /* for all atoms in the dihedral */
1517 GMX_ASSERT(vsnral
!= 0, "If we've seen a vsite before, we know how many constructing atoms it had");
1518 GMX_ASSERT(first_atoms
!= NULL
, "If we've seen a vsite before, we know what its first atom index was");
1519 atom
= ps
->param
[i
].a
[k
];
1520 if (vsite_type
[atom
] == NOTSET
&& vsite_type
[atom
] != F_VSITEN
)
1522 /* vsnral will be set here, we don't get here with nvsite==0 */
1524 for (m
= 0; (m
< vsnral
) && !bUsed
; m
++)
1526 if (atom
== first_atoms
[m
])
1536 fprintf(debug
, "unused atom in dih: %d\n", atom
+1);
1544 ps
->param
[kept_i
] = ps
->param
[i
];
1549 if (ps
->nr
!= kept_i
)
1551 fprintf(stderr
, "Removed %4d %15ss with virtual sites, %5d left\n",
1552 ps
->nr
-kept_i
, interaction_function
[cftype
].longname
, kept_i
);
1557 void clean_vsite_bondeds(t_params
*plist
, int natoms
, gmx_bool bRmVSiteBds
)
1559 int i
, k
, nvsite
, ftype
, vsite
, parnr
;
1562 at2vsitecon_t
*at2vc
;
1564 pindex
= nullptr; /* avoid warnings */
1565 /* make vsite_type array */
1566 snew(vsite_type
, natoms
);
1567 for (i
= 0; i
< natoms
; i
++)
1569 vsite_type
[i
] = NOTSET
;
1572 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1574 if (interaction_function
[ftype
].flags
& IF_VSITE
)
1576 nvsite
+= plist
[ftype
].nr
;
1578 while (i
< plist
[ftype
].nr
)
1580 vsite
= plist
[ftype
].param
[i
].ai();
1581 if (vsite_type
[vsite
] == NOTSET
)
1583 vsite_type
[vsite
] = ftype
;
1587 gmx_fatal(FARGS
, "multiple vsite constructions for atom %d", vsite
+1);
1589 if (ftype
== F_VSITEN
)
1591 while (i
< plist
[ftype
].nr
&& plist
[ftype
].param
[i
].ai() == vsite
)
1604 /* the rest only if we have virtual sites: */
1607 fprintf(stderr
, "Cleaning up constraints %swith virtual sites\n",
1608 bRmVSiteBds
? "and constant bonded interactions " : "");
1610 /* Make a reverse list to avoid ninteractions^2 operations */
1611 at2vc
= make_at2vsitecon(natoms
, plist
);
1613 snew(pindex
, natoms
);
1614 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1616 /* Here we skip VSITEN. In neary all practical use cases this
1617 * is not an issue, since VSITEN is intended for constructing
1618 * additional interaction sites, not for replacing normal atoms
1619 * with bonded interactions. Thus we do not expect constant
1620 * bonded interactions. If VSITEN does get used with constant
1621 * bonded interactions, these are not removed which only leads
1622 * to very minor extra computation and constant energy.
1623 * The only problematic case is onstraints between VSITEN
1624 * constructions with fixed distance (which is anyhow useless).
1625 * This will generate a fatal error in check_vsite_constraints.
1627 if ((interaction_function
[ftype
].flags
& IF_VSITE
) &&
1630 for (parnr
= 0; (parnr
< plist
[ftype
].nr
); parnr
++)
1632 k
= plist
[ftype
].param
[parnr
].ai();
1633 pindex
[k
].ftype
= ftype
;
1634 pindex
[k
].parnr
= parnr
;
1641 for (i
= 0; i
< natoms
; i
++)
1643 fprintf(debug
, "atom %d vsite_type %s\n", i
,
1644 vsite_type
[i
] == NOTSET
? "NOTSET" :
1645 interaction_function
[vsite_type
[i
]].name
);
1649 /* remove interactions that include virtual sites */
1650 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1652 if ( ( ( interaction_function
[ftype
].flags
& IF_BOND
) && bRmVSiteBds
) ||
1653 ( interaction_function
[ftype
].flags
& IF_CONSTRAINT
) )
1655 if (interaction_function
[ftype
].flags
& (IF_BTYPE
| IF_CONSTRAINT
) )
1657 clean_vsite_bonds (plist
, pindex
, ftype
, vsite_type
);
1659 else if (interaction_function
[ftype
].flags
& IF_ATYPE
)
1661 clean_vsite_angles(plist
, pindex
, ftype
, vsite_type
, at2vc
);
1663 else if ( (ftype
== F_PDIHS
) || (ftype
== F_IDIHS
) )
1665 clean_vsite_dihs (plist
, pindex
, ftype
, vsite_type
);
1669 /* check that no remaining constraints include virtual sites */
1670 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
1672 if (interaction_function
[ftype
].flags
& IF_CONSTRAINT
)
1674 check_vsite_constraints(plist
, ftype
, vsite_type
);
1678 done_at2vsitecon(natoms
, at2vc
);