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,2018 by the GROMACS development team.
7 * Copyright (c) 2019,2020,2021, by the GROMACS development team, led by
8 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
9 * and including many others, as listed in the AUTHORS file in the
10 * top-level source 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.
40 #include "vsite_parm.h"
49 #include "gromacs/gmxpreprocess/add_par.h"
50 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
51 #include "gromacs/gmxpreprocess/grompp_impl.h"
52 #include "gromacs/gmxpreprocess/notset.h"
53 #include "gromacs/gmxpreprocess/toputil.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/math/utilities.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/mdtypes/md_enums.h"
59 #include "gromacs/topology/ifunc.h"
60 #include "gromacs/topology/topology.h"
61 #include "gromacs/utility/arrayref.h"
62 #include "gromacs/utility/cstringutil.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/gmxassert.h"
65 #include "gromacs/utility/logger.h"
66 #include "gromacs/utility/smalloc.h"
67 #include "gromacs/utility/strconvert.h"
69 #include "hackblock.h"
73 * Data type to store information about bonded interactions for virtual sites.
75 class VsiteBondedInteraction
78 //! Constructor initializes datastructure.
79 VsiteBondedInteraction(gmx::ArrayRef<const int> atomIndex, real parameterValue) :
80 parameterValue_(parameterValue)
82 GMX_RELEASE_ASSERT(atomIndex.size() <= atomIndex_.size(),
83 "Cannot add more atom indices than maximum number");
84 std::array<int, 4>::iterator atomIndexIt = atomIndex_.begin();
85 for (const auto index : atomIndex)
87 *atomIndexIt++ = index;
91 //! Access the individual elements set for the vsite bonded interaction.
92 const int& ai() const { return atomIndex_[0]; }
93 const int& aj() const { return atomIndex_[1]; }
94 const int& ak() const { return atomIndex_[2]; }
95 const int& al() const { return atomIndex_[3]; }
97 const real& parameterValue() const { return parameterValue_; }
101 //! The distance value for this bonded interaction.
102 real parameterValue_;
103 //! Array of atom indices
104 std::array<int, 4> atomIndex_;
108 * Stores information about single virtual site bonded parameter.
110 struct VsiteBondParameter
112 //! Constructor initializes datastructure.
113 VsiteBondParameter(int ftype, const InteractionOfType& vsiteInteraction) :
114 ftype_(ftype), vsiteInteraction_(vsiteInteraction)
117 //! Function type for virtual site.
120 * Interaction type data.
122 * The datastructure should never be used in a case where the InteractionType
123 * used to construct it might go out of scope before this object, as it would cause
124 * the reference to type_ to dangle.
126 const InteractionOfType& vsiteInteraction_;
130 * Helper type for conversion of bonded parameters to virtual sites.
132 struct Atom2VsiteBond
134 //! Function type for conversion.
136 //! The vsite parameters in a list.
137 std::vector<VsiteBondParameter> vSiteBondedParameters;
140 //! Convenience type def for virtual site connections.
141 using Atom2VsiteConnection = std::vector<int>;
143 static int vsite_bond_nrcheck(int ftype)
147 if ((interaction_function[ftype].flags & (IF_BTYPE | IF_CONSTRAINT | IF_ATYPE)) || (ftype == F_IDIHS))
149 nrcheck = NRAL(ftype);
159 static void enter_bonded(int nratoms, std::vector<VsiteBondedInteraction>* bondeds, const InteractionOfType& type)
161 GMX_RELEASE_ASSERT(nratoms == type.atoms().ssize(), "Size of atom array must match");
162 bondeds->emplace_back(type.atoms(), type.c0());
166 * Wraps the datastructures for the different vsite bondeds.
168 struct AllVsiteBondedInteractions
171 std::vector<VsiteBondedInteraction> bonds;
173 std::vector<VsiteBondedInteraction> angles;
175 std::vector<VsiteBondedInteraction> dihedrals;
178 static AllVsiteBondedInteractions createVsiteBondedInformation(int nrat,
179 gmx::ArrayRef<const int> atoms,
180 gmx::ArrayRef<const Atom2VsiteBond> at2vb)
182 AllVsiteBondedInteractions allVsiteBondeds;
183 for (int k = 0; k < nrat; k++)
185 for (const auto& vsite : at2vb[atoms[k]].vSiteBondedParameters)
187 int ftype = vsite.ftype_;
188 const InteractionOfType& type = vsite.vsiteInteraction_;
189 int nrcheck = vsite_bond_nrcheck(ftype);
190 /* abuse nrcheck to see if we're adding bond, angle or idih */
193 case 2: enter_bonded(nrcheck, &allVsiteBondeds.bonds, type); break;
194 case 3: enter_bonded(nrcheck, &allVsiteBondeds.angles, type); break;
195 case 4: enter_bonded(nrcheck, &allVsiteBondeds.dihedrals, type); break;
199 return allVsiteBondeds;
202 static std::vector<Atom2VsiteBond> make_at2vsitebond(int natoms, gmx::ArrayRef<InteractionsOfType> plist)
206 std::vector<Atom2VsiteBond> at2vb(natoms);
209 for (int ftype = 0; (ftype < F_NRE); ftype++)
211 if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN)
213 for (int i = 0; (i < gmx::ssize(plist[ftype])); i++)
215 gmx::ArrayRef<const int> atoms = plist[ftype].interactionTypes[i].atoms();
216 for (int j = 0; j < NRAL(ftype); j++)
218 bVSI[atoms[j]] = TRUE;
224 for (int ftype = 0; (ftype < F_NRE); ftype++)
226 int nrcheck = vsite_bond_nrcheck(ftype);
229 for (int i = 0; (i < gmx::ssize(plist[ftype])); i++)
231 gmx::ArrayRef<const int> aa = plist[ftype].interactionTypes[i].atoms();
232 for (int j = 0; j < nrcheck; j++)
236 at2vb[aa[j]].vSiteBondedParameters.emplace_back(
237 ftype, plist[ftype].interactionTypes[i]);
248 static std::vector<Atom2VsiteConnection> make_at2vsitecon(int natoms, gmx::ArrayRef<InteractionsOfType> plist)
250 std::vector<bool> bVSI(natoms);
251 std::vector<Atom2VsiteConnection> at2vc(natoms);
253 for (int ftype = 0; (ftype < F_NRE); ftype++)
255 if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN)
257 for (int i = 0; (i < gmx::ssize(plist[ftype])); i++)
259 gmx::ArrayRef<const int> atoms = plist[ftype].interactionTypes[i].atoms();
260 for (int j = 0; j < NRAL(ftype); j++)
262 bVSI[atoms[j]] = TRUE;
268 for (int ftype = 0; (ftype < F_NRE); ftype++)
270 if (interaction_function[ftype].flags & IF_CONSTRAINT)
272 for (int i = 0; (i < gmx::ssize(plist[ftype])); i++)
274 int ai = plist[ftype].interactionTypes[i].ai();
275 int aj = plist[ftype].interactionTypes[i].aj();
276 if (bVSI[ai] && bVSI[aj])
278 /* Store forward direction */
279 at2vc[ai].emplace_back(aj);
280 /* Store backward direction */
281 at2vc[aj].emplace_back(ai);
290 static void print_bad(FILE* fp,
291 gmx::ArrayRef<const VsiteBondedInteraction> bonds,
292 gmx::ArrayRef<const VsiteBondedInteraction> angles,
293 gmx::ArrayRef<const VsiteBondedInteraction> idihs)
297 fprintf(fp, "bonds:");
298 for (const auto& bond : bonds)
300 fprintf(fp, " %d-%d (%g)", bond.ai() + 1, bond.aj() + 1, bond.parameterValue());
306 fprintf(fp, "angles:");
307 for (const auto& angle : angles)
309 fprintf(fp, " %d-%d-%d (%g)", angle.ai() + 1, angle.aj() + 1, angle.ak() + 1, angle.parameterValue());
315 fprintf(fp, "idihs:");
316 for (const auto& idih : idihs)
324 idih.parameterValue());
330 static void printInteractionOfType(FILE* fp, int ftype, int i, const InteractionOfType& type)
333 static int prev_ftype = NOTSET;
334 static int prev_i = NOTSET;
336 if ((ftype != prev_ftype) || (i != prev_i))
342 fprintf(fp, "(%d) plist[%s].param[%d]", pass, interaction_function[ftype].name, i);
343 gmx::ArrayRef<const real> forceParam = type.forceParam();
344 for (int j = 0; j < NRFP(ftype); j++)
346 fprintf(fp, ".c[%d]=%g ", j, forceParam[j]);
352 static real get_bond_length(gmx::ArrayRef<const VsiteBondedInteraction> bonds, t_iatom ai, t_iatom aj)
357 for (const auto& bond : bonds)
359 /* check both ways */
360 if (((ai == bond.ai()) && (aj == bond.aj())) || ((ai == bond.aj()) && (aj == bond.ai())))
362 bondlen = bond.parameterValue(); /* note: parameterValue might be NOTSET */
369 static real get_angle(gmx::ArrayRef<const VsiteBondedInteraction> angles, t_iatom ai, t_iatom aj, t_iatom ak)
374 for (const auto& ang : angles)
376 /* check both ways */
377 if (((ai == ang.ai()) && (aj == ang.aj()) && (ak == ang.ak()))
378 || ((ai == ang.ak()) && (aj == ang.aj()) && (ak == ang.ai())))
380 angle = gmx::c_deg2Rad * ang.parameterValue();
387 static const char* get_atomtype_name_AB(t_atom* atom, PreprocessingAtomTypes* atypes)
389 const char* name = *atypes->atomNameFromAtomType(atom->type);
391 /* When using the decoupling option, atom types are changed
392 * to decoupled for the non-bonded interactions, but the virtual
393 * sites constructions should be based on the "normal" interactions.
394 * So we return the state B atom type name if the state A atom
395 * type is the decoupled one. We should actually check for the atom
396 * type number, but that's not passed here. So we check for
397 * the decoupled atom type name. This should not cause trouble
398 * as this code is only used for topologies with v-sites without
399 * parameters generated by pdb2gmx.
401 if (strcmp(name, "decoupled") == 0)
403 name = *atypes->atomNameFromAtomType(atom->typeB);
409 static bool calc_vsite3_param(PreprocessingAtomTypes* atypes,
410 InteractionOfType* vsite,
412 gmx::ArrayRef<const VsiteBondedInteraction> bonds,
413 gmx::ArrayRef<const VsiteBondedInteraction> angles)
415 /* i = virtual site | ,k
416 * j = 1st bonded heavy atom | i-j
417 * k,l = 2nd bonded atoms | `l
421 real bjk, bjl, a = -1, b = -1;
422 /* check if this is part of a NH3 , NH2-umbrella or CH3 group,
423 * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
424 bXH3 = ((gmx::equalCaseInsensitive(get_atomtype_name_AB(&at->atom[vsite->ak()], atypes), "MNH", 3))
425 && (gmx::equalCaseInsensitive(
426 get_atomtype_name_AB(&at->atom[vsite->al()], atypes), "MNH", 3)))
427 || ((gmx::equalCaseInsensitive(
428 get_atomtype_name_AB(&at->atom[vsite->ak()], atypes), "MCH3", 4))
429 && (gmx::equalCaseInsensitive(
430 get_atomtype_name_AB(&at->atom[vsite->al()], atypes), "MCH3", 4)));
432 bjk = get_bond_length(bonds, vsite->aj(), vsite->ak());
433 bjl = get_bond_length(bonds, vsite->aj(), vsite->al());
434 bError = (bjk == NOTSET) || (bjl == NOTSET);
437 /* now we get some XH2/XH3 group specific construction */
438 /* note: we call the heavy atom 'C' and the X atom 'N' */
439 real bMM, bCM, bCN, bNH, aCNH, dH, rH, dM, rM;
442 /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
443 bError = bError || (bjk != bjl);
445 /* the X atom (C or N) in the XH2/XH3 group is the first after the masses: */
446 aN = std::max(vsite->ak(), vsite->al()) + 1;
448 /* get common bonds */
449 bMM = get_bond_length(bonds, vsite->ak(), vsite->al());
451 bCN = get_bond_length(bonds, vsite->aj(), aN);
452 bError = bError || (bMM == NOTSET) || (bCN == NOTSET);
454 /* calculate common things */
456 dM = std::sqrt(gmx::square(bCM) - gmx::square(rM));
458 /* are we dealing with the X atom? */
459 if (vsite->ai() == aN)
461 /* this is trivial */
462 a = b = 0.5 * bCN / dM;
466 /* get other bondlengths and angles: */
467 bNH = get_bond_length(bonds, aN, vsite->ai());
468 aCNH = get_angle(angles, vsite->aj(), aN, vsite->ai());
469 bError = bError || (bNH == NOTSET) || (aCNH == NOTSET);
472 dH = bCN - bNH * std::cos(aCNH);
473 rH = bNH * std::sin(aCNH);
475 a = 0.5 * (dH / dM + rH / rM);
476 b = 0.5 * (dH / dM - rH / rM);
482 "calc_vsite3_param not implemented for the general case "
486 vsite->setForceParameter(0, a);
487 vsite->setForceParameter(1, b);
492 static bool calc_vsite3fd_param(InteractionOfType* vsite,
493 gmx::ArrayRef<const VsiteBondedInteraction> bonds,
494 gmx::ArrayRef<const VsiteBondedInteraction> angles)
496 /* i = virtual site | ,k
497 * j = 1st bonded heavy atom | i-j
498 * k,l = 2nd bonded atoms | `l
502 real bij, bjk, bjl, aijk, aijl, rk, rl;
504 bij = get_bond_length(bonds, vsite->ai(), vsite->aj());
505 bjk = get_bond_length(bonds, vsite->aj(), vsite->ak());
506 bjl = get_bond_length(bonds, vsite->aj(), vsite->al());
507 aijk = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
508 aijl = get_angle(angles, vsite->ai(), vsite->aj(), vsite->al());
509 bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) || (aijk == NOTSET) || (aijl == NOTSET);
511 rk = bjk * std::sin(aijk);
512 rl = bjl * std::sin(aijl);
513 vsite->setForceParameter(0, rk / (rk + rl));
514 vsite->setForceParameter(1, -bij);
519 static bool calc_vsite3fad_param(InteractionOfType* vsite,
520 gmx::ArrayRef<const VsiteBondedInteraction> bonds,
521 gmx::ArrayRef<const VsiteBondedInteraction> angles)
523 /* i = virtual site |
524 * j = 1st bonded heavy atom | i-j
525 * k = 2nd bonded heavy atom | `k-l
526 * l = 3d bonded heavy atom |
529 bool bSwapParity, bError;
532 bSwapParity = (vsite->c1() == -1);
534 bij = get_bond_length(bonds, vsite->ai(), vsite->aj());
535 aijk = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
536 bError = (bij == NOTSET) || (aijk == NOTSET);
538 vsite->setForceParameter(1, bij);
539 vsite->setForceParameter(0, gmx::c_rad2Deg * aijk);
543 vsite->setForceParameter(0, 360 - vsite->c0());
549 static bool calc_vsite3out_param(PreprocessingAtomTypes* atypes,
550 InteractionOfType* vsite,
552 gmx::ArrayRef<const VsiteBondedInteraction> bonds,
553 gmx::ArrayRef<const VsiteBondedInteraction> angles)
555 /* i = virtual site | ,k
556 * j = 1st bonded heavy atom | i-j
557 * k,l = 2nd bonded atoms | `l
558 * NOTE: i is out of the j-k-l plane!
561 bool bXH3, bError, bSwapParity;
562 real bij, bjk, bjl, aijk, aijl, akjl, pijk, pijl, a, b, c;
564 /* check if this is part of a NH2-umbrella, NH3 or CH3 group,
565 * i.e. if atom k and l are dummy masses (MNH* or MCH3*) */
566 bXH3 = ((gmx::equalCaseInsensitive(get_atomtype_name_AB(&at->atom[vsite->ak()], atypes), "MNH", 3))
567 && (gmx::equalCaseInsensitive(
568 get_atomtype_name_AB(&at->atom[vsite->al()], atypes), "MNH", 3)))
569 || ((gmx::equalCaseInsensitive(
570 get_atomtype_name_AB(&at->atom[vsite->ak()], atypes), "MCH3", 4))
571 && (gmx::equalCaseInsensitive(
572 get_atomtype_name_AB(&at->atom[vsite->al()], atypes), "MCH3", 4)));
574 /* check if construction parity must be swapped */
575 bSwapParity = (vsite->c1() == -1);
577 bjk = get_bond_length(bonds, vsite->aj(), vsite->ak());
578 bjl = get_bond_length(bonds, vsite->aj(), vsite->al());
579 bError = (bjk == NOTSET) || (bjl == NOTSET);
582 /* now we get some XH3 group specific construction */
583 /* note: we call the heavy atom 'C' and the X atom 'N' */
584 real bMM, bCM, bCN, bNH, aCNH, dH, rH, rHx, rHy, dM, rM;
587 /* check if bonds from heavy atom (j) to dummy masses (k,l) are equal: */
588 bError = bError || (bjk != bjl);
590 /* the X atom (C or N) in the XH3 group is the first after the masses: */
591 aN = std::max(vsite->ak(), vsite->al()) + 1;
593 /* get all bondlengths and angles: */
594 bMM = get_bond_length(bonds, vsite->ak(), vsite->al());
596 bCN = get_bond_length(bonds, vsite->aj(), aN);
597 bNH = get_bond_length(bonds, aN, vsite->ai());
598 aCNH = get_angle(angles, vsite->aj(), aN, vsite->ai());
599 bError = bError || (bMM == NOTSET) || (bCN == NOTSET) || (bNH == NOTSET) || (aCNH == NOTSET);
602 dH = bCN - bNH * std::cos(aCNH);
603 rH = bNH * std::sin(aCNH);
604 /* we assume the H's are symmetrically distributed */
605 rHx = rH * std::cos(gmx::c_deg2Rad * 30);
606 rHy = rH * std::sin(gmx::c_deg2Rad * 30);
608 dM = std::sqrt(gmx::square(bCM) - gmx::square(rM));
609 a = 0.5 * ((dH / dM) - (rHy / rM));
610 b = 0.5 * ((dH / dM) + (rHy / rM));
611 c = rHx / (2 * dM * rM);
615 /* this is the general construction */
617 bij = get_bond_length(bonds, vsite->ai(), vsite->aj());
618 aijk = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
619 aijl = get_angle(angles, vsite->ai(), vsite->aj(), vsite->al());
620 akjl = get_angle(angles, vsite->ak(), vsite->aj(), vsite->al());
621 bError = bError || (bij == NOTSET) || (aijk == NOTSET) || (aijl == NOTSET) || (akjl == NOTSET);
623 pijk = std::cos(aijk) * bij;
624 pijl = std::cos(aijl) * bij;
625 a = (pijk + (pijk * std::cos(akjl) - pijl) * std::cos(akjl) / gmx::square(std::sin(akjl))) / bjk;
626 b = (pijl + (pijl * std::cos(akjl) - pijk) * std::cos(akjl) / gmx::square(std::sin(akjl))) / bjl;
627 c = -std::sqrt(gmx::square(bij)
628 - (gmx::square(pijk) - 2 * pijk * pijl * std::cos(akjl) + gmx::square(pijl))
629 / gmx::square(std::sin(akjl)))
630 / (bjk * bjl * std::sin(akjl));
633 vsite->setForceParameter(0, a);
634 vsite->setForceParameter(1, b);
637 vsite->setForceParameter(2, -c);
641 vsite->setForceParameter(2, c);
646 static bool calc_vsite4fd_param(InteractionOfType* vsite,
647 gmx::ArrayRef<const VsiteBondedInteraction> bonds,
648 gmx::ArrayRef<const VsiteBondedInteraction> angles,
649 const gmx::MDLogger& logger)
651 /* i = virtual site | ,k
652 * j = 1st bonded heavy atom | i-j-m
653 * k,l,m = 2nd bonded atoms | `l
657 real bij, bjk, bjl, bjm, aijk, aijl, aijm, akjm, akjl;
658 real pk, pl, pm, cosakl, cosakm, sinakl, sinakm, cl, cm;
660 bij = get_bond_length(bonds, vsite->ai(), vsite->aj());
661 bjk = get_bond_length(bonds, vsite->aj(), vsite->ak());
662 bjl = get_bond_length(bonds, vsite->aj(), vsite->al());
663 bjm = get_bond_length(bonds, vsite->aj(), vsite->am());
664 aijk = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
665 aijl = get_angle(angles, vsite->ai(), vsite->aj(), vsite->al());
666 aijm = get_angle(angles, vsite->ai(), vsite->aj(), vsite->am());
667 akjm = get_angle(angles, vsite->ak(), vsite->aj(), vsite->am());
668 akjl = get_angle(angles, vsite->ak(), vsite->aj(), vsite->al());
669 bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) || (bjm == NOTSET) || (aijk == NOTSET)
670 || (aijl == NOTSET) || (aijm == NOTSET) || (akjm == NOTSET) || (akjl == NOTSET);
674 pk = bjk * std::sin(aijk);
675 pl = bjl * std::sin(aijl);
676 pm = bjm * std::sin(aijm);
677 cosakl = (std::cos(akjl) - std::cos(aijk) * std::cos(aijl)) / (std::sin(aijk) * std::sin(aijl));
678 cosakm = (std::cos(akjm) - std::cos(aijk) * std::cos(aijm)) / (std::sin(aijk) * std::sin(aijm));
679 if (cosakl < -1 || cosakl > 1 || cosakm < -1 || cosakm > 1)
681 GMX_LOG(logger.warning)
683 .appendTextFormatted(
684 "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f",
686 gmx::c_rad2Deg * aijk,
687 gmx::c_rad2Deg * aijl,
688 gmx::c_rad2Deg * aijm);
690 "invalid construction in calc_vsite4fd for atom %d: "
691 "cosakl=%f, cosakm=%f\n",
696 sinakl = std::sqrt(1 - gmx::square(cosakl));
697 sinakm = std::sqrt(1 - gmx::square(cosakm));
699 /* note: there is a '+' because of the way the sines are calculated */
700 cl = -pk / (pl * cosakl - pk + pl * sinakl * (pm * cosakm - pk) / (pm * sinakm));
701 cm = -pk / (pm * cosakm - pk + pm * sinakm * (pl * cosakl - pk) / (pl * sinakl));
703 vsite->setForceParameter(0, cl);
704 vsite->setForceParameter(1, cm);
705 vsite->setForceParameter(2, -bij);
712 static bool calc_vsite4fdn_param(InteractionOfType* vsite,
713 gmx::ArrayRef<const VsiteBondedInteraction> bonds,
714 gmx::ArrayRef<const VsiteBondedInteraction> angles,
715 const gmx::MDLogger& logger)
717 /* i = virtual site | ,k
718 * j = 1st bonded heavy atom | i-j-m
719 * k,l,m = 2nd bonded atoms | `l
723 real bij, bjk, bjl, bjm, aijk, aijl, aijm;
724 real pk, pl, pm, a, b;
726 bij = get_bond_length(bonds, vsite->ai(), vsite->aj());
727 bjk = get_bond_length(bonds, vsite->aj(), vsite->ak());
728 bjl = get_bond_length(bonds, vsite->aj(), vsite->al());
729 bjm = get_bond_length(bonds, vsite->aj(), vsite->am());
730 aijk = get_angle(angles, vsite->ai(), vsite->aj(), vsite->ak());
731 aijl = get_angle(angles, vsite->ai(), vsite->aj(), vsite->al());
732 aijm = get_angle(angles, vsite->ai(), vsite->aj(), vsite->am());
734 bError = (bij == NOTSET) || (bjk == NOTSET) || (bjl == NOTSET) || (bjm == NOTSET)
735 || (aijk == NOTSET) || (aijl == NOTSET) || (aijm == NOTSET);
740 /* Calculate component of bond j-k along the direction i-j */
741 pk = -bjk * std::cos(aijk);
743 /* Calculate component of bond j-l along the direction i-j */
744 pl = -bjl * std::cos(aijl);
746 /* Calculate component of bond j-m along the direction i-j */
747 pm = -bjm * std::cos(aijm);
749 if (fabs(pl) < 1000 * GMX_REAL_MIN || fabs(pm) < 1000 * GMX_REAL_MIN)
751 GMX_LOG(logger.warning)
753 .appendTextFormatted(
754 "virtual site %d: angle ijk = %f, angle ijl = %f, angle ijm = %f",
756 gmx::c_rad2Deg * aijk,
757 gmx::c_rad2Deg * aijl,
758 gmx::c_rad2Deg * aijm);
760 "invalid construction in calc_vsite4fdn for atom %d: "
770 vsite->setForceParameter(0, a);
771 vsite->setForceParameter(1, b);
772 vsite->setForceParameter(2, bij);
779 int set_vsites(bool bVerbose,
781 PreprocessingAtomTypes* atypes,
782 gmx::ArrayRef<InteractionsOfType> plist,
783 const gmx::MDLogger& logger)
792 /* Make a reverse list to avoid ninteractions^2 operations */
793 std::vector<Atom2VsiteBond> at2vb = make_at2vsitebond(atoms->nr, plist);
795 for (ftype = 0; (ftype < F_NRE); ftype++)
797 if (interaction_function[ftype].flags & IF_VSITE)
799 nvsite += plist[ftype].size();
801 if (ftype == F_VSITEN)
803 /* We don't calculate parameters for VSITEN */
809 for (auto& param : plist[ftype].interactionTypes)
811 /* check if all parameters are set */
813 gmx::ArrayRef<const real> forceParam = param.forceParam();
814 for (int j = 0; (j < NRFP(ftype)) && bSet; j++)
816 bSet = forceParam[j] != NOTSET;
821 fprintf(debug, "bSet=%s ", gmx::boolToString(bSet));
822 printInteractionOfType(debug, ftype, i, plist[ftype].interactionTypes[i]);
826 if (bVerbose && bFirst)
830 .appendTextFormatted("Calculating parameters for virtual sites");
835 /* now set the vsite parameters: */
836 AllVsiteBondedInteractions allVsiteBondeds =
837 createVsiteBondedInformation(NRAL(ftype), param.atoms(), at2vb);
841 "Found %zu bonds, %zu angles and %zu idihs "
842 "for virtual site %d (%s)\n",
843 allVsiteBondeds.bonds.size(),
844 allVsiteBondeds.angles.size(),
845 allVsiteBondeds.dihedrals.size(),
847 interaction_function[ftype].longname);
849 allVsiteBondeds.bonds,
850 allVsiteBondeds.angles,
851 allVsiteBondeds.dihedrals);
856 bERROR = calc_vsite3_param(
857 atypes, ¶m, atoms, allVsiteBondeds.bonds, allVsiteBondeds.angles);
860 bERROR = calc_vsite3fd_param(
861 ¶m, allVsiteBondeds.bonds, allVsiteBondeds.angles);
864 bERROR = calc_vsite3fad_param(
865 ¶m, allVsiteBondeds.bonds, allVsiteBondeds.angles);
868 bERROR = calc_vsite3out_param(
869 atypes, ¶m, atoms, allVsiteBondeds.bonds, allVsiteBondeds.angles);
872 bERROR = calc_vsite4fd_param(
873 ¶m, allVsiteBondeds.bonds, allVsiteBondeds.angles, logger);
876 bERROR = calc_vsite4fdn_param(
877 ¶m, allVsiteBondeds.bonds, allVsiteBondeds.angles, logger);
881 "Automatic parameter generation not supported "
883 interaction_function[ftype].longname,
890 "Automatic parameter generation not supported "
891 "for %s atom %d for this bonding configuration",
892 interaction_function[ftype].longname,
903 void set_vsites_ptype(bool bVerbose, gmx_moltype_t* molt, const gmx::MDLogger& logger)
911 .appendTextFormatted("Setting particle type to V for virtual sites");
913 for (ftype = 0; ftype < F_NRE; ftype++)
915 InteractionList* il = &molt->ilist[ftype];
916 if (interaction_function[ftype].flags & IF_VSITE)
918 const int nra = interaction_function[ftype].nratoms;
919 const int nrd = il->size();
920 gmx::ArrayRef<const int> ia = il->iatoms;
926 .appendTextFormatted("doing %d %s virtual sites",
928 interaction_function[ftype].longname);
931 for (i = 0; (i < nrd);)
933 /* The virtual site */
934 int avsite = ia[i + 1];
935 molt->atoms.atom[avsite].ptype = ParticleType::VSite;
944 * Convenience typedef for linking function type to parameter numbers.
946 * The entries in this datastructure are valid if the particle participates in
947 * a virtual site interaction and has a valid vsite function type other than VSITEN.
948 * \todo Change to remove empty constructor when gmx::compat::optional is available.
950 class VsiteAtomMapping
953 //! Only construct with all information in place or nothing
954 VsiteAtomMapping(int functionType, int interactionIndex) :
955 functionType_(functionType), interactionIndex_(interactionIndex)
958 VsiteAtomMapping() : functionType_(-1), interactionIndex_(-1) {}
959 //! Get function type.
960 const int& functionType() const { return functionType_; }
961 //! Get parameter number.
962 const int& interactionIndex() const { return interactionIndex_; }
965 //! Function type for the linked parameter.
967 //! The linked parameter.
968 int interactionIndex_;
971 static void check_vsite_constraints(gmx::ArrayRef<InteractionsOfType> plist,
973 const int vsite_type[],
974 const gmx::MDLogger& logger)
977 for (const auto& param : plist[cftype].interactionTypes)
979 gmx::ArrayRef<const int> atoms = param.atoms();
980 for (int k = 0; k < 2; k++)
983 if (vsite_type[atom] != NOTSET)
987 .appendTextFormatted(
988 "ERROR: Cannot have constraint (%d-%d) with virtual site (%d)",
998 gmx_fatal(FARGS, "There were %d virtual sites involved in constraints", n);
1002 static void clean_vsite_bonds(gmx::ArrayRef<InteractionsOfType> plist,
1003 gmx::ArrayRef<const VsiteAtomMapping> pindex,
1005 const int vsite_type[],
1006 const gmx::MDLogger& logger)
1009 int nconverted, nremoved;
1010 int oatom, at1, at2;
1011 bool bKeep, bRemove, bAllFD;
1012 InteractionsOfType* ps;
1014 if (cftype == F_CONNBONDS)
1019 ps = &(plist[cftype]);
1023 for (auto bond = ps->interactionTypes.begin(); bond != ps->interactionTypes.end();)
1026 const int* first_atoms = nullptr;
1031 /* check if all virtual sites are constructed from the same atoms */
1033 gmx::ArrayRef<const int> atoms = bond->atoms();
1034 for (int k = 0; (k < 2) && !bKeep && !bRemove; k++)
1036 /* for all atoms in the bond */
1037 int atom = atoms[k];
1038 if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN)
1041 bool bThisFD = ((pindex[atom].functionType() == F_VSITE3FD)
1042 || (pindex[atom].functionType() == F_VSITE3FAD)
1043 || (pindex[atom].functionType() == F_VSITE4FD)
1044 || (pindex[atom].functionType() == F_VSITE4FDN));
1045 bool bThisOUT = ((pindex[atom].functionType() == F_VSITE3OUT)
1046 && ((interaction_function[cftype].flags & IF_CONSTRAINT) != 0U));
1047 bAllFD = bAllFD && bThisFD;
1048 if (bThisFD || bThisOUT)
1050 oatom = atoms[1 - k]; /* the other atom */
1051 if (vsite_type[oatom] == NOTSET
1053 == plist[pindex[atom].functionType()]
1054 .interactionTypes[pindex[atom].interactionIndex()]
1057 /* if the other atom isn't a vsite, and it is AI */
1067 /* TODO This fragment, and corresponding logic in
1068 clean_vsite_angles and clean_vsite_dihs should
1069 be refactored into a common function */
1072 /* if this is the first vsite we encounter then
1073 store construction atoms */
1074 /* TODO This would be nicer to implement with
1075 a C++ "vector view" class" with an
1076 STL-container-like interface. */
1077 vsnral = NRAL(pindex[atom].functionType()) - 1;
1078 first_atoms = plist[pindex[atom].functionType()]
1079 .interactionTypes[pindex[atom].interactionIndex()]
1086 GMX_ASSERT(vsnral != 0, "nvsite > 1 must have vsnral != 0");
1087 GMX_ASSERT(first_atoms != nullptr,
1088 "nvsite > 1 must have first_atoms != NULL");
1089 /* if it is not the first then
1090 check if this vsite is constructed from the same atoms */
1091 if (vsnral == NRAL(pindex[atom].functionType()) - 1)
1093 for (int m = 0; (m < vsnral) && !bKeep; m++)
1097 bool bPresent = false;
1098 atoms = plist[pindex[atom].functionType()]
1099 .interactionTypes[pindex[atom].interactionIndex()]
1103 for (int n = 0; (n < vsnral) && !bPresent; n++)
1105 if (atoms[m] == first_atoms[n])
1131 /* if we have no virtual sites in this bond, keep it */
1137 /* TODO This loop and the corresponding loop in
1138 check_vsite_angles should be refactored into a common
1140 /* check if all non-vsite atoms are used in construction: */
1141 bool bFirstTwo = true;
1142 for (int k = 0; (k < 2) && !bKeep; k++) /* for all atoms in the bond */
1144 int atom = atoms[k];
1145 if (vsite_type[atom] == NOTSET)
1148 for (int m = 0; (m < vsnral) && !bUsed; m++)
1150 GMX_ASSERT(first_atoms != nullptr,
1151 "If we've seen a vsite before, we know what its first atom "
1154 if (atom == first_atoms[m])
1157 bFirstTwo = bFirstTwo && m < 2;
1167 if (!(bAllFD && bFirstTwo))
1169 /* Two atom bonded interactions include constraints.
1170 * We need to remove constraints between vsite pairs that have
1171 * a fixed distance due to being constructed from the same
1172 * atoms, since this can be numerically unstable.
1174 for (int m = 0; m < vsnral && !bKeep; m++) /* all constr. atoms */
1176 at1 = first_atoms[m];
1177 at2 = first_atoms[(m + 1) % vsnral];
1178 bool bPresent = false;
1179 for (ftype = 0; ftype < F_NRE; ftype++)
1181 if (interaction_function[ftype].flags & IF_CONSTRAINT)
1183 for (auto entry = plist[ftype].interactionTypes.begin();
1184 (entry != plist[ftype].interactionTypes.end()) && !bPresent;
1187 /* all constraints until one matches */
1188 bPresent = (((entry->ai() == at1) && (entry->aj() == at2))
1189 || ((entry->ai() == at2) && (entry->aj() == at1)));
1205 else if (IS_CHEMBOND(cftype))
1207 plist[F_CONNBONDS].interactionTypes.emplace_back(*bond);
1208 bond = ps->interactionTypes.erase(bond);
1213 bond = ps->interactionTypes.erase(bond);
1220 GMX_LOG(logger.info)
1222 .appendTextFormatted("Removed %4d %15ss with virtual sites, %zu left",
1224 interaction_function[cftype].longname,
1229 GMX_LOG(logger.info)
1231 .appendTextFormatted(
1232 "Converted %4d %15ss with virtual sites to connections, %zu left",
1234 interaction_function[cftype].longname,
1239 GMX_LOG(logger.info)
1241 .appendTextFormatted(
1242 "Warning: removed %d %ss with vsite with %s construction\n"
1243 " This vsite construction does not guarantee constant "
1245 " If the constructions were generated by pdb2gmx ignore "
1248 interaction_function[cftype].longname,
1249 interaction_function[F_VSITE3OUT].longname);
1253 static void clean_vsite_angles(gmx::ArrayRef<InteractionsOfType> plist,
1254 gmx::ArrayRef<VsiteAtomMapping> pindex,
1256 const int vsite_type[],
1257 gmx::ArrayRef<const Atom2VsiteConnection> at2vc,
1258 const gmx::MDLogger& logger)
1261 InteractionsOfType* ps;
1263 ps = &(plist[cftype]);
1264 int oldSize = ps->size();
1265 for (auto angle = ps->interactionTypes.begin(); angle != ps->interactionTypes.end();)
1268 const int* first_atoms = nullptr;
1271 bool bAll3FAD = true;
1272 /* check if all virtual sites are constructed from the same atoms */
1274 gmx::ArrayRef<const int> atoms = angle->atoms();
1275 for (int k = 0; (k < 3) && !bKeep; k++) /* for all atoms in the angle */
1277 int atom = atoms[k];
1278 if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN)
1281 bAll3FAD = bAll3FAD && (pindex[atom].functionType() == F_VSITE3FAD);
1284 /* store construction atoms of first vsite */
1285 vsnral = NRAL(pindex[atom].functionType()) - 1;
1286 first_atoms = plist[pindex[atom].functionType()]
1287 .interactionTypes[pindex[atom].interactionIndex()]
1294 GMX_ASSERT(vsnral != 0,
1295 "If we've seen a vsite before, we know how many constructing atoms "
1298 first_atoms != nullptr,
1299 "If we've seen a vsite before, we know what its first atom index was");
1300 /* check if this vsite is constructed from the same atoms */
1301 if (vsnral == NRAL(pindex[atom].functionType()) - 1)
1303 for (int m = 0; (m < vsnral) && !bKeep; m++)
1305 const int* subAtoms;
1307 bool bPresent = false;
1308 subAtoms = plist[pindex[atom].functionType()]
1309 .interactionTypes[pindex[atom].interactionIndex()]
1313 for (int n = 0; (n < vsnral) && !bPresent; n++)
1315 if (subAtoms[m] == first_atoms[n])
1334 /* keep all angles with no virtual sites in them or
1335 with virtual sites with more than 3 constr. atoms */
1336 if (nvsite == 0 && vsnral > 3)
1341 /* check if all non-vsite atoms are used in construction: */
1342 bool bFirstTwo = true;
1343 for (int k = 0; (k < 3) && !bKeep; k++) /* for all atoms in the angle */
1346 if (vsite_type[atom] == NOTSET)
1349 for (int m = 0; (m < vsnral) && !bUsed; m++)
1352 first_atoms != nullptr,
1353 "If we've seen a vsite before, we know what its first atom index was");
1355 if (atom == first_atoms[m])
1358 bFirstTwo = bFirstTwo && m < 2;
1368 if (!(bAll3FAD && bFirstTwo))
1370 /* check if all constructing atoms are constrained together */
1371 for (int m = 0; m < vsnral && !bKeep; m++) /* all constr. atoms */
1373 at1 = first_atoms[m];
1374 at2 = first_atoms[(m + 1) % vsnral];
1375 bool bPresent = false;
1376 auto found = std::find(at2vc[at1].begin(), at2vc[at1].end(), at2);
1377 if (found != at2vc[at1].end())
1394 angle = ps->interactionTypes.erase(angle);
1398 if (oldSize != gmx::ssize(*ps))
1400 GMX_LOG(logger.info)
1402 .appendTextFormatted("Removed %4zu %15ss with virtual sites, %zu left",
1403 oldSize - ps->size(),
1404 interaction_function[cftype].longname,
1409 static void clean_vsite_dihs(gmx::ArrayRef<InteractionsOfType> plist,
1410 gmx::ArrayRef<const VsiteAtomMapping> pindex,
1412 const int vsite_type[],
1413 const gmx::MDLogger& logger)
1415 InteractionsOfType* ps;
1417 ps = &(plist[cftype]);
1419 int oldSize = ps->size();
1420 for (auto dih = ps->interactionTypes.begin(); dih != ps->interactionTypes.end();)
1423 const int* first_atoms = nullptr;
1426 gmx::ArrayRef<const int> atoms = dih->atoms();
1428 /* check if all virtual sites are constructed from the same atoms */
1430 for (int k = 0; (k < 4) && !bKeep; k++) /* for all atoms in the dihedral */
1433 if (vsite_type[atom] != NOTSET && vsite_type[atom] != F_VSITEN)
1437 /* store construction atoms of first vsite */
1438 vsnral = NRAL(pindex[atom].functionType()) - 1;
1439 first_atoms = plist[pindex[atom].functionType()]
1440 .interactionTypes[pindex[atom].interactionIndex()]
1447 GMX_ASSERT(vsnral != 0,
1448 "If we've seen a vsite before, we know how many constructing atoms "
1451 first_atoms != nullptr,
1452 "If we've seen a vsite before, we know what its first atom index was");
1453 /* check if this vsite is constructed from the same atoms */
1454 if (vsnral == NRAL(pindex[atom].functionType()) - 1)
1456 for (int m = 0; (m < vsnral) && !bKeep; m++)
1458 const int* subAtoms;
1460 bool bPresent = false;
1461 subAtoms = plist[pindex[atom].functionType()]
1462 .interactionTypes[pindex[atom].interactionIndex()]
1466 for (int n = 0; (n < vsnral) && !bPresent; n++)
1468 if (subAtoms[m] == first_atoms[n])
1480 /* TODO clean_site_bonds and _angles do this increment
1481 at the top of the loop. Refactor this for
1487 /* keep all dihedrals with no virtual sites in them */
1493 /* check if all atoms in dihedral are either virtual sites, or used in
1494 construction of virtual sites. If so, keep it, if not throw away: */
1495 for (int k = 0; (k < 4) && !bKeep; k++) /* for all atoms in the dihedral */
1497 GMX_ASSERT(vsnral != 0,
1498 "If we've seen a vsite before, we know how many constructing atoms it had");
1499 GMX_ASSERT(first_atoms != nullptr,
1500 "If we've seen a vsite before, we know what its first atom index was");
1502 if (vsite_type[atom] == NOTSET)
1504 /* vsnral will be set here, we don't get here with nvsite==0 */
1506 for (int m = 0; (m < vsnral) && !bUsed; m++)
1508 if (atom == first_atoms[m])
1526 dih = ps->interactionTypes.erase(dih);
1530 if (oldSize != gmx::ssize(*ps))
1532 GMX_LOG(logger.info)
1534 .appendTextFormatted("Removed %4zu %15ss with virtual sites, %zu left",
1535 oldSize - ps->size(),
1536 interaction_function[cftype].longname,
1541 // TODO use gmx::compat::optional for pindex.
1542 void clean_vsite_bondeds(gmx::ArrayRef<InteractionsOfType> plist,
1545 const gmx::MDLogger& logger)
1549 std::vector<VsiteAtomMapping> pindex;
1550 std::vector<Atom2VsiteConnection> at2vc;
1552 /* make vsite_type array */
1553 snew(vsite_type, natoms);
1554 for (int i = 0; i < natoms; i++)
1556 vsite_type[i] = NOTSET;
1559 for (int ftype = 0; ftype < F_NRE; ftype++)
1561 if (interaction_function[ftype].flags & IF_VSITE)
1563 nvsite += plist[ftype].size();
1565 while (i < gmx::ssize(plist[ftype]))
1567 vsite = plist[ftype].interactionTypes[i].ai();
1568 if (vsite_type[vsite] == NOTSET)
1570 vsite_type[vsite] = ftype;
1574 gmx_fatal(FARGS, "multiple vsite constructions for atom %d", vsite + 1);
1576 if (ftype == F_VSITEN)
1578 while (i < gmx::ssize(plist[ftype]) && plist[ftype].interactionTypes[i].ai() == vsite)
1591 /* the rest only if we have virtual sites: */
1594 GMX_LOG(logger.info)
1596 .appendTextFormatted("Cleaning up constraints %swith virtual sites",
1597 bRmVSiteBds ? "and constant bonded interactions " : "");
1599 /* Make a reverse list to avoid ninteractions^2 operations */
1600 at2vc = make_at2vsitecon(natoms, plist);
1602 pindex.resize(natoms);
1603 for (int ftype = 0; ftype < F_NRE; ftype++)
1605 /* Here we skip VSITEN. In neary all practical use cases this
1606 * is not an issue, since VSITEN is intended for constructing
1607 * additional interaction sites, not for replacing normal atoms
1608 * with bonded interactions. Thus we do not expect constant
1609 * bonded interactions. If VSITEN does get used with constant
1610 * bonded interactions, these are not removed which only leads
1611 * to very minor extra computation and constant energy.
1612 * The only problematic case is onstraints between VSITEN
1613 * constructions with fixed distance (which is anyhow useless).
1614 * This will generate a fatal error in check_vsite_constraints.
1616 if ((interaction_function[ftype].flags & IF_VSITE) && ftype != F_VSITEN)
1618 for (gmx::index interactionIndex = 0; interactionIndex < gmx::ssize(plist[ftype]);
1621 int k = plist[ftype].interactionTypes[interactionIndex].ai();
1622 pindex[k] = VsiteAtomMapping(ftype, interactionIndex);
1627 /* remove interactions that include virtual sites */
1628 for (int ftype = 0; ftype < F_NRE; ftype++)
1630 if (((interaction_function[ftype].flags & IF_BOND) && bRmVSiteBds)
1631 || (interaction_function[ftype].flags & IF_CONSTRAINT))
1633 if (interaction_function[ftype].flags & (IF_BTYPE | IF_CONSTRAINT))
1635 clean_vsite_bonds(plist, pindex, ftype, vsite_type, logger);
1637 else if (interaction_function[ftype].flags & IF_ATYPE)
1639 clean_vsite_angles(plist, pindex, ftype, vsite_type, at2vc, logger);
1641 else if ((ftype == F_PDIHS) || (ftype == F_IDIHS))
1643 clean_vsite_dihs(plist, pindex, ftype, vsite_type, logger);
1647 /* check that no remaining constraints include virtual sites */
1648 for (int ftype = 0; ftype < F_NRE; ftype++)
1650 if (interaction_function[ftype].flags & IF_CONSTRAINT)
1652 check_vsite_constraints(plist, ftype, vsite_type, logger);