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,2016,2017 by the GROMACS development team.
7 * Copyright (c) 2018,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 "gen_vsite.h"
52 #include "gromacs/fileio/pdbio.h"
53 #include "gromacs/gmxpreprocess/add_par.h"
54 #include "gromacs/gmxpreprocess/fflibutil.h"
55 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
56 #include "gromacs/gmxpreprocess/grompp_impl.h"
57 #include "gromacs/gmxpreprocess/notset.h"
58 #include "gromacs/gmxpreprocess/toputil.h"
59 #include "gromacs/math/functions.h"
60 #include "gromacs/math/units.h"
61 #include "gromacs/math/vec.h"
62 #include "gromacs/mdtypes/md_enums.h"
63 #include "gromacs/topology/ifunc.h"
64 #include "gromacs/topology/residuetypes.h"
65 #include "gromacs/topology/symtab.h"
66 #include "gromacs/utility/basedefinitions.h"
67 #include "gromacs/utility/cstringutil.h"
68 #include "gromacs/utility/fatalerror.h"
69 #include "gromacs/utility/futil.h"
70 #include "gromacs/utility/real.h"
71 #include "gromacs/utility/smalloc.h"
73 #include "hackblock.h"
77 #define OPENDIR '[' /* starting sign for directive */
78 #define CLOSEDIR ']' /* ending sign for directive */
80 /*! \libinternal \brief
81 * The configuration describing a virtual site.
83 struct VirtualSiteConfiguration
86 * Explicit constructor.
88 * \param[in] type Atomtype for vsite configuration.
89 * \param[in] planar Is the input conf planar.
90 * \param[in] nhyd How many hydrogens are in the configuration.
91 * \param[in] nextheavy Type of bonded heavy atom.
92 * \param[in] dummy What kind of dummy is used in the vsite.
94 explicit VirtualSiteConfiguration(const std::string& type,
97 const std::string& nextheavy,
98 const std::string& dummy) :
102 nextHeavyType(nextheavy),
106 //! Type for the XH3/XH2 atom.
107 std::string atomtype;
108 /*! \brief Is the configuration planar?
110 * If true, the atomtype above and the three connected
111 * ones are in a planar geometry. The two next entries
112 * are undefined in that case.
114 bool isplanar = false;
115 //! cnumber of connected hydrogens.
117 //! Type for the heavy atom bonded to XH2/XH3.
118 std::string nextHeavyType;
119 //! The type of MNH* or MCH3* dummy mass to use.
120 std::string dummyMass;
124 /*!\libinternal \brief
125 * Virtual site topology datastructure.
127 * Structure to represent average bond and angles values in vsite aromatic
128 * residues. Note that these are NOT necessarily the bonds and angles from the
129 * forcefield; many forcefields (like Amber, OPLS) have some inherent strain in
130 * 5-rings (i.e. the sum of angles is !=540, but impropers keep it planar)
132 struct VirtualSiteTopology
135 * Explicit constructor
137 * \param[in] name Residue name.
139 explicit VirtualSiteTopology(const std::string& name) : resname(name) {}
142 //! Helper struct for single bond in virtual site.
143 struct VirtualSiteBond
146 * Explicit constructor
148 * \param[in] a1 First atom name.
149 * \param[in] a2 Second atom name.
150 * \param[in] v Value for distance.
152 VirtualSiteBond(const std::string& a1, const std::string& a2, real v) :
162 //! Distance value between atoms.
165 //! Container of all bonds in virtual site.
166 std::vector<VirtualSiteBond> bond;
167 //! Helper struct for single angle in virtual site.
168 struct VirtualSiteAngle
171 * Explicit constructor
173 * \param[in] a1 First atom name.
174 * \param[in] a2 Second atom name.
175 * \param[in] a3 Third atom name.
176 * \param[in] v Value for angle.
178 VirtualSiteAngle(const std::string& a1, const std::string& a2, const std::string& a3, real v) :
194 //! Container for all angles in virtual site.
195 std::vector<VirtualSiteAngle> angle;
213 typedef char t_dirname[STRLEN];
215 static const t_dirname ddb_dirnames[DDB_DIR_NR] = { "CH3", "NH3", "NH2", "PHE", "TYR",
216 "TRP", "HISA", "HISB", "HISH" };
218 static int ddb_name2dir(char* name)
220 /* Translate a directive name to the number of the directive.
221 * HID/HIE/HIP names are translated to the ones we use in Gromacs.
228 for (i = 0; i < DDB_DIR_NR && index < 0; i++)
230 if (!gmx_strcasecmp(name, ddb_dirnames[i]))
240 static void read_vsite_database(const char* ddbname,
241 std::vector<VirtualSiteConfiguration>* vsiteconflist,
242 std::vector<VirtualSiteTopology>* vsitetoplist)
244 /* This routine is a quick hack to fix the problem with hardcoded atomtypes
245 * and aromatic vsite parameters by reading them from a ff???.vsd file.
247 * The file can contain sections [ NH3 ], [ CH3 ], [ NH2 ], and ring residue names.
248 * For the NH3 and CH3 section each line has three fields. The first is the atomtype
249 * (nb: not bonded type) of the N/C atom to be replaced, the second field is
250 * the type of the next heavy atom it is bonded to, and the third field the type
251 * of dummy mass that will be used for this group.
253 * If the NH2 group planar (sp2 N) a different vsite construct is used, so in this
254 * case the second field should just be the word planar.
261 char s1[MAXNAME], s2[MAXNAME], s3[MAXNAME], s4[MAXNAME];
263 gmx::FilePtr ddb = gmx::openLibraryFile(ddbname);
267 while (fgets2(pline, STRLEN - 2, ddb.get()) != nullptr)
269 strip_comment(pline);
271 if (strlen(pline) > 0)
273 if (pline[0] == OPENDIR)
275 strncpy(dirstr, pline + 1, STRLEN - 2);
276 if ((ch = strchr(dirstr, CLOSEDIR)) != nullptr)
282 if (!gmx_strcasecmp(dirstr, "HID") || !gmx_strcasecmp(dirstr, "HISD"))
284 sprintf(dirstr, "HISA");
286 else if (!gmx_strcasecmp(dirstr, "HIE") || !gmx_strcasecmp(dirstr, "HISE"))
288 sprintf(dirstr, "HISB");
290 else if (!gmx_strcasecmp(dirstr, "HIP"))
292 sprintf(dirstr, "HISH");
295 curdir = ddb_name2dir(dirstr);
298 gmx_fatal(FARGS, "Invalid directive %s in vsite database %s", dirstr, ddbname);
306 gmx_fatal(FARGS, "First entry in vsite database must be a directive.\n");
311 int numberOfSites = sscanf(pline, "%s%s%s", s1, s2, s3);
312 std::string s1String = s1;
313 std::string s2String = s2;
314 std::string s3String = s3;
315 if (numberOfSites < 3 && gmx::equalCaseInsensitive(s2String, "planar"))
317 VirtualSiteConfiguration newVsiteConf(s1String, true, 2, "0", "0");
318 vsiteconflist->push_back(newVsiteConf);
320 else if (numberOfSites == 3)
322 VirtualSiteConfiguration newVsiteConf(s1String, false, -1, s2String, s3String);
323 if (curdir == DDB_NH2)
325 newVsiteConf.nHydrogens = 2;
329 newVsiteConf.nHydrogens = 3;
331 vsiteconflist->push_back(newVsiteConf);
335 gmx_fatal(FARGS, "Not enough directives in vsite database line: %s\n", pline);
346 const auto found = std::find_if(
347 vsitetoplist->begin(), vsitetoplist->end(), [&dirstr](const auto& entry) {
348 return gmx::equalCaseInsensitive(dirstr, entry.resname);
350 /* Allocate a new topology entry if this is a new residue */
351 if (found == vsitetoplist->end())
353 vsitetoplist->push_back(VirtualSiteTopology(dirstr));
355 int numberOfSites = sscanf(pline, "%s%s%s%s", s1, s2, s3, s4);
356 std::string s1String = s1;
357 std::string s2String = s2;
358 std::string s3String = s3;
360 if (numberOfSites == 3)
363 vsitetoplist->back().bond.emplace_back(s1String, s2String, strtod(s3, nullptr));
365 else if (numberOfSites == 4)
368 vsitetoplist->back().angle.emplace_back(s1String, s2String, s3String,
369 strtod(s4, nullptr));
375 "Need 3 or 4 values to specify bond/angle values in %s: %s\n",
382 "Didn't find a case for directive %s in read_vsite_database\n", dirstr);
389 static int nitrogen_is_planar(gmx::ArrayRef<const VirtualSiteConfiguration> vsiteconflist,
390 const std::string& atomtype)
392 /* Return 1 if atomtype exists in database list and is planar, 0 if not,
393 * and -1 if not found.
397 std::find_if(vsiteconflist.begin(), vsiteconflist.end(), [&atomtype](const auto& entry) {
398 return (gmx::equalCaseInsensitive(entry.atomtype, atomtype) && entry.nHydrogens == 2);
400 if (found != vsiteconflist.end())
402 res = static_cast<int>(found->isplanar);
412 static std::string get_dummymass_name(gmx::ArrayRef<const VirtualSiteConfiguration> vsiteconflist,
413 const std::string& atom,
414 const std::string& nextheavy)
416 /* Return the dummy mass name if found, or NULL if not set in ddb database */
417 const auto found = std::find_if(
418 vsiteconflist.begin(), vsiteconflist.end(), [&atom, &nextheavy](const auto& entry) {
419 return (gmx::equalCaseInsensitive(atom, entry.atomtype)
420 && gmx::equalCaseInsensitive(nextheavy, entry.nextHeavyType));
422 if (found != vsiteconflist.end())
424 return found->dummyMass;
433 static real get_ddb_bond(gmx::ArrayRef<const VirtualSiteTopology> vsitetop,
434 const std::string& res,
435 const std::string& atom1,
436 const std::string& atom2)
438 const auto found = std::find_if(vsitetop.begin(), vsitetop.end(), [&res](const auto& entry) {
439 return gmx::equalCaseInsensitive(res, entry.resname);
442 if (found == vsitetop.end())
444 gmx_fatal(FARGS, "No vsite information for residue %s found in vsite database.\n", res.c_str());
446 const auto foundBond =
447 std::find_if(found->bond.begin(), found->bond.end(), [&atom1, &atom2](const auto& entry) {
448 return ((atom1 == entry.atom1 && atom2 == entry.atom2)
449 || (atom1 == entry.atom2 && atom2 == entry.atom1));
451 if (foundBond == found->bond.end())
453 gmx_fatal(FARGS, "Couldnt find bond %s-%s for residue %s in vsite database.\n",
454 atom1.c_str(), atom2.c_str(), res.c_str());
457 return foundBond->value;
461 static real get_ddb_angle(gmx::ArrayRef<const VirtualSiteTopology> vsitetop,
462 const std::string& res,
463 const std::string& atom1,
464 const std::string& atom2,
465 const std::string& atom3)
467 const auto found = std::find_if(vsitetop.begin(), vsitetop.end(), [&res](const auto& entry) {
468 return gmx::equalCaseInsensitive(res, entry.resname);
471 if (found == vsitetop.end())
473 gmx_fatal(FARGS, "No vsite information for residue %s found in vsite database.\n", res.c_str());
475 const auto foundAngle = std::find_if(
476 found->angle.begin(), found->angle.end(), [&atom1, &atom2, &atom3](const auto& entry) {
477 return ((atom1 == entry.atom1 && atom2 == entry.atom2 && atom3 == entry.atom3)
478 || (atom1 == entry.atom3 && atom2 == entry.atom2 && atom3 == entry.atom1)
479 || (atom1 == entry.atom2 && atom2 == entry.atom1 && atom3 == entry.atom3)
480 || (atom1 == entry.atom3 && atom2 == entry.atom1 && atom3 == entry.atom2));
483 if (foundAngle == found->angle.end())
485 gmx_fatal(FARGS, "Couldnt find angle %s-%s-%s for residue %s in vsite database.\n",
486 atom1.c_str(), atom2.c_str(), atom3.c_str(), res.c_str());
489 return foundAngle->value;
493 static void count_bonds(int atom,
494 InteractionsOfType* psb,
503 int heavy, other, nrb, nrH, nrhv;
505 /* find heavy atom bound to this hydrogen */
507 for (auto parm = psb->interactionTypes.begin();
508 (parm != psb->interactionTypes.end()) && (heavy == NOTSET); parm++)
510 if (parm->ai() == atom)
514 else if (parm->aj() == atom)
521 gmx_fatal(FARGS, "unbound hydrogen atom %d", atom + 1);
523 /* find all atoms bound to heavy atom */
528 for (const auto& parm : psb->interactionTypes)
530 if (parm.ai() == heavy)
534 else if (parm.aj() == heavy)
541 if (is_hydrogen(*(atomname[other])))
548 heavies[nrhv] = other;
561 print_bonds(FILE* fp, int o2n[], int nrHatoms, const int Hatoms[], int Heavy, int nrheavies, const int heavies[])
565 fprintf(fp, "Found: %d Hatoms: ", nrHatoms);
566 for (i = 0; i < nrHatoms; i++)
568 fprintf(fp, " %d", o2n[Hatoms[i]] + 1);
570 fprintf(fp, "; %d Heavy atoms: %d", nrheavies + 1, o2n[Heavy] + 1);
571 for (i = 0; i < nrheavies; i++)
573 fprintf(fp, " %d", o2n[heavies[i]] + 1);
578 static int get_atype(int atom, t_atoms* at, gmx::ArrayRef<const PreprocessResidue> rtpFFDB, ResidueType* rt)
583 if (at->atom[atom].m != 0.0F)
585 type = at->atom[atom].type;
589 /* get type from rtpFFDB */
590 auto localPpResidue = getDatabaseEntry(*(at->resinfo[at->atom[atom].resind].name), rtpFFDB);
591 bNterm = rt->namedResidueHasType(*(at->resinfo[at->atom[atom].resind].name), "Protein")
592 && (at->atom[atom].resind == 0);
593 int j = search_jtype(*localPpResidue, *(at->atomname[atom]), bNterm);
594 type = localPpResidue->atom[j].type;
599 static int vsite_nm2type(const char* name, PreprocessingAtomTypes* atype)
603 tp = atype->atomTypeFromName(name);
606 gmx_fatal(FARGS, "Dummy mass type (%s) not found in atom type database", name);
612 static real get_amass(int atom, t_atoms* at, gmx::ArrayRef<const PreprocessResidue> rtpFFDB, ResidueType* rt)
617 if (at->atom[atom].m != 0.0F)
619 mass = at->atom[atom].m;
623 /* get mass from rtpFFDB */
624 auto localPpResidue = getDatabaseEntry(*(at->resinfo[at->atom[atom].resind].name), rtpFFDB);
625 bNterm = rt->namedResidueHasType(*(at->resinfo[at->atom[atom].resind].name), "Protein")
626 && (at->atom[atom].resind == 0);
627 int j = search_jtype(*localPpResidue, *(at->atomname[atom]), bNterm);
628 mass = localPpResidue->atom[j].m;
633 static void my_add_param(InteractionsOfType* plist, int ai, int aj, real b)
635 static real c[MAXFORCEPARAM] = { NOTSET, NOTSET, NOTSET, NOTSET, NOTSET, NOTSET };
638 add_param(plist, ai, aj, c, nullptr);
641 static void add_vsites(gmx::ArrayRef<InteractionsOfType> plist,
649 int other, moreheavy;
651 for (int i = 0; i < nrHatoms; i++)
653 int ftype = vsite_type[Hatoms[i]];
654 /* Errors in setting the vsite_type should really be caugth earlier,
655 * because here it's not possible to print any useful error message.
656 * But it's still better to print a message than to segfault.
660 gmx_incons("Undetected error in setting up virtual sites");
662 bool bSwapParity = (ftype < 0);
663 vsite_type[Hatoms[i]] = ftype = abs(ftype);
664 if (ftype == F_BONDS)
666 if ((nrheavies != 1) && (nrHatoms != 1))
669 "cannot make constraint in add_vsites for %d heavy "
670 "atoms and %d hydrogen atoms",
671 nrheavies, nrHatoms);
673 my_add_param(&(plist[F_CONSTRNC]), Hatoms[i], heavies[0], NOTSET);
684 gmx_fatal(FARGS, "Not enough heavy atoms (%d) for %s (min 3)",
685 nrheavies + 1, interaction_function[vsite_type[Hatoms[i]]].name);
687 add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], heavies[1], bSwapParity);
693 moreheavy = heavies[1];
697 /* find more heavy atoms */
698 other = moreheavy = NOTSET;
699 for (auto parm = plist[F_BONDS].interactionTypes.begin();
700 (parm != plist[F_BONDS].interactionTypes.end()) && (moreheavy == NOTSET);
703 if (parm->ai() == heavies[0])
707 else if (parm->aj() == heavies[0])
711 if ((other != NOTSET) && (other != Heavy))
716 if (moreheavy == NOTSET)
718 gmx_fatal(FARGS, "Unbound molecule part %d-%d", Heavy + 1, Hatoms[0] + 1);
721 add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], moreheavy, bSwapParity);
728 gmx_fatal(FARGS, "Not enough heavy atoms (%d) for %s (min 4)",
729 nrheavies + 1, interaction_function[vsite_type[Hatoms[i]]].name);
731 add_vsite4_atoms(&plist[ftype], Hatoms[0], Heavy, heavies[0], heavies[1], heavies[2]);
735 gmx_fatal(FARGS, "can't use add_vsites for interaction function %s",
736 interaction_function[vsite_type[Hatoms[i]]].name);
742 #define ANGLE_6RING (DEG2RAD * 120)
744 /* cosine rule: a^2 = b^2 + c^2 - 2 b c cos(alpha) */
745 /* get a^2 when a, b and alpha are given: */
746 #define cosrule(b, c, alpha) (gmx::square(b) + gmx::square(c) - 2 * (b) * (c)*std::cos(alpha))
747 /* get cos(alpha) when a, b and c are given: */
748 #define acosrule(a, b, c) ((gmx::square(b) + gmx::square(c) - gmx::square(a)) / (2 * (b) * (c)))
750 static int gen_vsites_6ring(t_atoms* at,
752 gmx::ArrayRef<InteractionsOfType> plist,
760 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
778 real a, b, dCGCE, tmp1, tmp2, mtot, mG, mrest;
780 /* CG, CE1 and CE2 stay and each get a part of the total mass,
781 * so the c-o-m stays the same.
788 gmx_incons("Generating vsites on 6-rings");
792 /* constraints between CG, CE1 and CE2: */
793 dCGCE = std::sqrt(cosrule(bond_cc, bond_cc, ANGLE_6RING));
794 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE);
795 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE2], dCGCE);
796 my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atCE2], dCGCE);
798 /* rest will be vsite3 */
801 for (i = 0; i < (bDoZ ? atNR : atHZ); i++)
803 mtot += at->atom[ats[i]].m;
804 if (i != atCG && i != atCE1 && i != atCE2 && (bDoZ || (i != atHZ && i != atCZ)))
806 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
807 (*vsite_type)[ats[i]] = F_VSITE3;
811 /* Distribute mass so center-of-mass stays the same.
812 * The center-of-mass in the call is defined with x=0 at
813 * the CE1-CE2 bond and y=0 at the line from CG to the middle of CE1-CE2 bond.
815 xCG = -bond_cc + bond_cc * std::cos(ANGLE_6RING);
817 mG = at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = xcom * mtot / xCG;
819 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB = at->atom[ats[atCE2]].m =
820 at->atom[ats[atCE2]].mB = mrest / 2;
822 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
823 tmp1 = dCGCE * std::sin(ANGLE_6RING * 0.5);
824 tmp2 = bond_cc * std::cos(0.5 * ANGLE_6RING) + tmp1;
826 a = b = -bond_ch / tmp1;
828 add_vsite3_param(&plist[F_VSITE3], ats[atHE1], ats[atCE1], ats[atCE2], ats[atCG], a, b);
829 add_vsite3_param(&plist[F_VSITE3], ats[atHE2], ats[atCE2], ats[atCE1], ats[atCG], a, b);
830 /* CD1, CD2 and CZ: */
832 add_vsite3_param(&plist[F_VSITE3], ats[atCD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
833 add_vsite3_param(&plist[F_VSITE3], ats[atCD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
836 add_vsite3_param(&plist[F_VSITE3], ats[atCZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
838 /* HD1, HD2 and HZ: */
839 a = b = (bond_ch + tmp2) / tmp1;
840 add_vsite3_param(&plist[F_VSITE3], ats[atHD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
841 add_vsite3_param(&plist[F_VSITE3], ats[atHD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
844 add_vsite3_param(&plist[F_VSITE3], ats[atHZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
850 static int gen_vsites_phe(t_atoms* at,
852 gmx::ArrayRef<InteractionsOfType> plist,
855 gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
857 real bond_cc, bond_ch;
860 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
877 /* Aromatic rings have 6-fold symmetry, so we only need one bond length.
878 * (angle is always 120 degrees).
880 bond_cc = get_ddb_bond(vsitetop, "PHE", "CD1", "CE1");
881 bond_ch = get_ddb_bond(vsitetop, "PHE", "CD1", "HD1");
883 x[atCG] = -bond_cc + bond_cc * std::cos(ANGLE_6RING);
885 x[atHD1] = x[atCD1] + bond_ch * std::cos(ANGLE_6RING);
887 x[atHE1] = x[atCE1] - bond_ch * std::cos(ANGLE_6RING);
892 x[atCZ] = bond_cc * std::cos(0.5 * ANGLE_6RING);
893 x[atHZ] = x[atCZ] + bond_ch;
896 for (i = 0; i < atNR; i++)
898 xcom += x[i] * at->atom[ats[i]].m;
899 mtot += at->atom[ats[i]].m;
903 return gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, TRUE);
907 calc_vsite3_param(real xd, real yd, real xi, real yi, real xj, real yj, real xk, real yk, real* a, real* b)
909 /* determine parameters by solving the equation system, since we know the
910 * virtual site coordinates here.
912 real dx_ij, dx_ik, dy_ij, dy_ik;
919 *a = ((xd - xi) * dy_ik - dx_ik * (yd - yi)) / (dx_ij * dy_ik - dx_ik * dy_ij);
920 *b = (yd - yi - (*a) * dy_ij) / dy_ik;
924 static int gen_vsites_trp(PreprocessingAtomTypes* atype,
925 std::vector<gmx::RVec>* newx,
927 char*** newatomname[],
929 int* newvsite_type[],
933 gmx::ArrayRef<const gmx::RVec> x,
937 gmx::ArrayRef<InteractionsOfType> plist,
941 gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
944 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
965 /* weights for determining the COM's of both rings (M1 and M2): */
966 real mw[NMASS][atNR] = { { 0, 1, 1, 1, 0.5, 1, 1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0 },
967 { 0, 0, 0, 0, 0.5, 0, 0, 0.5, 1, 1, 1, 1, 1, 1, 1, 1 } };
969 real xi[atNR], yi[atNR];
970 real xcom[NMASS], ycom[NMASS], alpha;
971 real b_CD2_CE2, b_NE1_CE2, b_CG_CD2, b_CH2_HH2, b_CE2_CZ2;
972 real b_NE1_HE1, b_CD2_CE3, b_CE3_CZ3, b_CB_CG;
973 real b_CZ2_CH2, b_CZ2_HZ2, b_CD1_HD1, b_CE3_HE3;
974 real b_CG_CD1, b_CZ3_HZ3;
975 real a_NE1_CE2_CD2, a_CE2_CD2_CG, a_CB_CG_CD2, a_CE2_CD2_CE3;
976 real a_CD2_CG_CD1, a_CE2_CZ2_HZ2, a_CZ2_CH2_HH2;
977 real a_CD2_CE2_CZ2, a_CD2_CE3_CZ3, a_CE3_CZ3_HZ3, a_CG_CD1_HD1;
978 real a_CE2_CZ2_CH2, a_HE1_NE1_CE2, a_CD2_CE3_HE3;
979 int atM[NMASS], tpM, i, i0, j, nvsite;
980 real mM[NMASS], dCBM1, dCBM2, dM1M2;
982 rvec r_ij, r_ik, t1, t2;
987 gmx_incons("atom types in gen_vsites_trp");
989 /* Get geometry from database */
990 b_CD2_CE2 = get_ddb_bond(vsitetop, "TRP", "CD2", "CE2");
991 b_NE1_CE2 = get_ddb_bond(vsitetop, "TRP", "NE1", "CE2");
992 b_CG_CD1 = get_ddb_bond(vsitetop, "TRP", "CG", "CD1");
993 b_CG_CD2 = get_ddb_bond(vsitetop, "TRP", "CG", "CD2");
994 b_CB_CG = get_ddb_bond(vsitetop, "TRP", "CB", "CG");
995 b_CE2_CZ2 = get_ddb_bond(vsitetop, "TRP", "CE2", "CZ2");
996 b_CD2_CE3 = get_ddb_bond(vsitetop, "TRP", "CD2", "CE3");
997 b_CE3_CZ3 = get_ddb_bond(vsitetop, "TRP", "CE3", "CZ3");
998 b_CZ2_CH2 = get_ddb_bond(vsitetop, "TRP", "CZ2", "CH2");
1000 b_CD1_HD1 = get_ddb_bond(vsitetop, "TRP", "CD1", "HD1");
1001 b_CZ2_HZ2 = get_ddb_bond(vsitetop, "TRP", "CZ2", "HZ2");
1002 b_NE1_HE1 = get_ddb_bond(vsitetop, "TRP", "NE1", "HE1");
1003 b_CH2_HH2 = get_ddb_bond(vsitetop, "TRP", "CH2", "HH2");
1004 b_CE3_HE3 = get_ddb_bond(vsitetop, "TRP", "CE3", "HE3");
1005 b_CZ3_HZ3 = get_ddb_bond(vsitetop, "TRP", "CZ3", "HZ3");
1007 a_NE1_CE2_CD2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "NE1", "CE2", "CD2");
1008 a_CE2_CD2_CG = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE2", "CD2", "CG");
1009 a_CB_CG_CD2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CB", "CG", "CD2");
1010 a_CD2_CG_CD1 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CD2", "CG", "CD1");
1012 a_CE2_CD2_CE3 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE2", "CD2", "CE3");
1013 a_CD2_CE2_CZ2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CD2", "CE2", "CZ2");
1014 a_CD2_CE3_CZ3 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CD2", "CE3", "CZ3");
1015 a_CE3_CZ3_HZ3 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE3", "CZ3", "HZ3");
1016 a_CZ2_CH2_HH2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CZ2", "CH2", "HH2");
1017 a_CE2_CZ2_HZ2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE2", "CZ2", "HZ2");
1018 a_CE2_CZ2_CH2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE2", "CZ2", "CH2");
1019 a_CG_CD1_HD1 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CG", "CD1", "HD1");
1020 a_HE1_NE1_CE2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "HE1", "NE1", "CE2");
1021 a_CD2_CE3_HE3 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CD2", "CE3", "HE3");
1023 /* Calculate local coordinates.
1024 * y-axis (x=0) is the bond CD2-CE2.
1025 * x-axis (y=0) is perpendicular to the bond CD2-CE2 and
1026 * intersects the middle of the bond.
1029 yi[atCD2] = -0.5 * b_CD2_CE2;
1032 yi[atCE2] = 0.5 * b_CD2_CE2;
1034 xi[atNE1] = -b_NE1_CE2 * std::sin(a_NE1_CE2_CD2);
1035 yi[atNE1] = yi[atCE2] - b_NE1_CE2 * std::cos(a_NE1_CE2_CD2);
1037 xi[atCG] = -b_CG_CD2 * std::sin(a_CE2_CD2_CG);
1038 yi[atCG] = yi[atCD2] + b_CG_CD2 * std::cos(a_CE2_CD2_CG);
1040 alpha = a_CE2_CD2_CG + M_PI - a_CB_CG_CD2;
1041 xi[atCB] = xi[atCG] - b_CB_CG * std::sin(alpha);
1042 yi[atCB] = yi[atCG] + b_CB_CG * std::cos(alpha);
1044 alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - M_PI;
1045 xi[atCD1] = xi[atCG] - b_CG_CD1 * std::sin(alpha);
1046 yi[atCD1] = yi[atCG] + b_CG_CD1 * std::cos(alpha);
1048 xi[atCE3] = b_CD2_CE3 * std::sin(a_CE2_CD2_CE3);
1049 yi[atCE3] = yi[atCD2] + b_CD2_CE3 * std::cos(a_CE2_CD2_CE3);
1051 xi[atCZ2] = b_CE2_CZ2 * std::sin(a_CD2_CE2_CZ2);
1052 yi[atCZ2] = yi[atCE2] - b_CE2_CZ2 * std::cos(a_CD2_CE2_CZ2);
1054 alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - M_PI;
1055 xi[atCZ3] = xi[atCE3] + b_CE3_CZ3 * std::sin(alpha);
1056 yi[atCZ3] = yi[atCE3] + b_CE3_CZ3 * std::cos(alpha);
1058 alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - M_PI;
1059 xi[atCH2] = xi[atCZ2] + b_CZ2_CH2 * std::sin(alpha);
1060 yi[atCH2] = yi[atCZ2] - b_CZ2_CH2 * std::cos(alpha);
1063 alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - a_CG_CD1_HD1;
1064 xi[atHD1] = xi[atCD1] - b_CD1_HD1 * std::sin(alpha);
1065 yi[atHD1] = yi[atCD1] + b_CD1_HD1 * std::cos(alpha);
1067 alpha = a_NE1_CE2_CD2 + M_PI - a_HE1_NE1_CE2;
1068 xi[atHE1] = xi[atNE1] - b_NE1_HE1 * std::sin(alpha);
1069 yi[atHE1] = yi[atNE1] - b_NE1_HE1 * std::cos(alpha);
1071 alpha = a_CE2_CD2_CE3 + M_PI - a_CD2_CE3_HE3;
1072 xi[atHE3] = xi[atCE3] + b_CE3_HE3 * std::sin(alpha);
1073 yi[atHE3] = yi[atCE3] + b_CE3_HE3 * std::cos(alpha);
1075 alpha = a_CD2_CE2_CZ2 + M_PI - a_CE2_CZ2_HZ2;
1076 xi[atHZ2] = xi[atCZ2] + b_CZ2_HZ2 * std::sin(alpha);
1077 yi[atHZ2] = yi[atCZ2] - b_CZ2_HZ2 * std::cos(alpha);
1079 alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - a_CZ2_CH2_HH2;
1080 xi[atHZ3] = xi[atCZ3] + b_CZ3_HZ3 * std::sin(alpha);
1081 yi[atHZ3] = yi[atCZ3] + b_CZ3_HZ3 * std::cos(alpha);
1083 alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - a_CE3_CZ3_HZ3;
1084 xi[atHH2] = xi[atCH2] + b_CH2_HH2 * std::sin(alpha);
1085 yi[atHH2] = yi[atCH2] - b_CH2_HH2 * std::cos(alpha);
1087 /* Calculate masses for each ring and put it on the dummy masses */
1088 for (j = 0; j < NMASS; j++)
1090 mM[j] = xcom[j] = ycom[j] = 0;
1092 for (i = 0; i < atNR; i++)
1096 for (j = 0; j < NMASS; j++)
1098 mM[j] += mw[j][i] * at->atom[ats[i]].m;
1099 xcom[j] += xi[i] * mw[j][i] * at->atom[ats[i]].m;
1100 ycom[j] += yi[i] * mw[j][i] * at->atom[ats[i]].m;
1104 for (j = 0; j < NMASS; j++)
1110 /* get dummy mass type */
1111 tpM = vsite_nm2type("MW", atype);
1112 /* make space for 2 masses: shift all atoms starting with CB */
1114 for (j = 0; j < NMASS; j++)
1116 atM[j] = i0 + *nadd + j;
1120 fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, (*o2n)[i0] + 1);
1123 for (j = i0; j < at->nr; j++)
1125 (*o2n)[j] = j + *nadd;
1127 newx->resize(at->nr + *nadd);
1128 srenew(*newatom, at->nr + *nadd);
1129 srenew(*newatomname, at->nr + *nadd);
1130 srenew(*newvsite_type, at->nr + *nadd);
1131 srenew(*newcgnr, at->nr + *nadd);
1132 for (j = 0; j < NMASS; j++)
1134 (*newatomname)[at->nr + *nadd - 1 - j] = nullptr;
1137 /* Dummy masses will be placed at the center-of-mass in each ring. */
1139 /* calc initial position for dummy masses in real (non-local) coordinates.
1140 * Cheat by using the routine to calculate virtual site parameters. It is
1141 * much easier when we have the coordinates expressed in terms of
1144 rvec_sub(x[ats[atCB]], x[ats[atCG]], r_ij);
1145 rvec_sub(x[ats[atCD2]], x[ats[atCG]], r_ik);
1146 calc_vsite3_param(xcom[0], ycom[0], xi[atCG], yi[atCG], xi[atCB], yi[atCB], xi[atCD2],
1150 rvec_add(t1, t2, t1);
1151 rvec_add(t1, x[ats[atCG]], (*newx)[atM[0]]);
1153 calc_vsite3_param(xcom[1], ycom[1], xi[atCG], yi[atCG], xi[atCB], yi[atCB], xi[atCD2],
1157 rvec_add(t1, t2, t1);
1158 rvec_add(t1, x[ats[atCG]], (*newx)[atM[1]]);
1160 /* set parameters for the masses */
1161 for (j = 0; j < NMASS; j++)
1163 sprintf(name, "MW%d", j + 1);
1164 (*newatomname)[atM[j]] = put_symtab(symtab, name);
1165 (*newatom)[atM[j]].m = (*newatom)[atM[j]].mB = mM[j];
1166 (*newatom)[atM[j]].q = (*newatom)[atM[j]].qB = 0.0;
1167 (*newatom)[atM[j]].type = (*newatom)[atM[j]].typeB = tpM;
1168 (*newatom)[atM[j]].ptype = eptAtom;
1169 (*newatom)[atM[j]].resind = at->atom[i0].resind;
1170 (*newatom)[atM[j]].elem[0] = 'M';
1171 (*newatom)[atM[j]].elem[1] = '\0';
1172 (*newvsite_type)[atM[j]] = NOTSET;
1173 (*newcgnr)[atM[j]] = (*cgnr)[i0];
1175 /* renumber cgnr: */
1176 for (i = i0; i < at->nr; i++)
1181 /* constraints between CB, M1 and M2 */
1182 /* 'add_shift' says which atoms won't be renumbered afterwards */
1183 dCBM1 = std::hypot(xcom[0] - xi[atCB], ycom[0] - yi[atCB]);
1184 dM1M2 = std::hypot(xcom[0] - xcom[1], ycom[0] - ycom[1]);
1185 dCBM2 = std::hypot(xcom[1] - xi[atCB], ycom[1] - yi[atCB]);
1186 my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift + atM[0], dCBM1);
1187 my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift + atM[1], dCBM2);
1188 my_add_param(&(plist[F_CONSTRNC]), add_shift + atM[0], add_shift + atM[1], dM1M2);
1190 /* rest will be vsite3 */
1192 for (i = 0; i < atNR; i++)
1196 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
1197 (*vsite_type)[ats[i]] = F_VSITE3;
1202 /* now define all vsites from M1, M2, CB, ie:
1203 r_d = r_M1 + a r_M1_M2 + b r_M1_CB */
1204 for (i = 0; i < atNR; i++)
1206 if ((*vsite_type)[ats[i]] == F_VSITE3)
1208 calc_vsite3_param(xi[i], yi[i], xcom[0], ycom[0], xcom[1], ycom[1], xi[atCB], yi[atCB],
1210 add_vsite3_param(&plist[F_VSITE3], ats[i], add_shift + atM[0], add_shift + atM[1],
1219 static int gen_vsites_tyr(PreprocessingAtomTypes* atype,
1220 std::vector<gmx::RVec>* newx,
1222 char*** newatomname[],
1224 int* newvsite_type[],
1228 gmx::ArrayRef<const gmx::RVec> x,
1232 gmx::ArrayRef<InteractionsOfType> plist,
1236 gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
1238 int nvsite, i, i0, j, atM, tpM;
1239 real dCGCE, dCEOH, dCGM, tmp1, a, b;
1240 real bond_cc, bond_ch, bond_co, bond_oh, angle_coh;
1246 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1264 /* CG, CE1, CE2 (as in general 6-ring) and OH and HH stay,
1265 rest gets virtualized.
1266 Now we have two linked triangles with one improper keeping them flat */
1267 if (atNR != nrfound)
1269 gmx_incons("Number of atom types in gen_vsites_tyr");
1272 /* Aromatic rings have 6-fold symmetry, so we only need one bond length
1273 * for the ring part (angle is always 120 degrees).
1275 bond_cc = get_ddb_bond(vsitetop, "TYR", "CD1", "CE1");
1276 bond_ch = get_ddb_bond(vsitetop, "TYR", "CD1", "HD1");
1277 bond_co = get_ddb_bond(vsitetop, "TYR", "CZ", "OH");
1278 bond_oh = get_ddb_bond(vsitetop, "TYR", "OH", "HH");
1279 angle_coh = DEG2RAD * get_ddb_angle(vsitetop, "TYR", "CZ", "OH", "HH");
1281 xi[atCG] = -bond_cc + bond_cc * std::cos(ANGLE_6RING);
1282 xi[atCD1] = -bond_cc;
1283 xi[atHD1] = xi[atCD1] + bond_ch * std::cos(ANGLE_6RING);
1285 xi[atHE1] = xi[atCE1] - bond_ch * std::cos(ANGLE_6RING);
1286 xi[atCD2] = xi[atCD1];
1287 xi[atHD2] = xi[atHD1];
1288 xi[atCE2] = xi[atCE1];
1289 xi[atHE2] = xi[atHE1];
1290 xi[atCZ] = bond_cc * std::cos(0.5 * ANGLE_6RING);
1291 xi[atOH] = xi[atCZ] + bond_co;
1294 for (i = 0; i < atOH; i++)
1296 xcom += xi[i] * at->atom[ats[i]].m;
1297 mtot += at->atom[ats[i]].m;
1301 /* first do 6 ring as default,
1302 except CZ (we'll do that different) and HZ (we don't have that): */
1303 nvsite = gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, FALSE);
1305 /* then construct CZ from the 2nd triangle */
1306 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
1307 a = b = 0.5 * bond_co / (bond_co - bond_cc * std::cos(ANGLE_6RING));
1308 add_vsite3_param(&plist[F_VSITE3], ats[atCZ], ats[atOH], ats[atCE1], ats[atCE2], a, b);
1309 at->atom[ats[atCZ]].m = at->atom[ats[atCZ]].mB = 0;
1311 /* constraints between CE1, CE2 and OH */
1312 dCGCE = std::sqrt(cosrule(bond_cc, bond_cc, ANGLE_6RING));
1313 dCEOH = std::sqrt(cosrule(bond_cc, bond_co, ANGLE_6RING));
1314 my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atOH], dCEOH);
1315 my_add_param(&(plist[F_CONSTRNC]), ats[atCE2], ats[atOH], dCEOH);
1317 /* We also want to constrain the angle C-O-H, but since CZ is constructed
1318 * we need to introduce a constraint to CG.
1319 * CG is much further away, so that will lead to instabilities in LINCS
1320 * when we constrain both CG-HH and OH-HH distances. Instead of requiring
1321 * the use of lincs_order=8 we introduce a dummy mass three times further
1322 * away from OH than HH. The mass is accordingly a third, with the remaining
1323 * 2/3 moved to OH. This shouldn't cause any problems since the forces will
1324 * apply to the HH constructed atom and not directly on the virtual mass.
1327 vdist = 2.0 * bond_oh;
1328 mM = at->atom[ats[atHH]].m / 2.0;
1329 at->atom[ats[atOH]].m += mM; /* add 1/2 of original H mass */
1330 at->atom[ats[atOH]].mB += mM; /* add 1/2 of original H mass */
1331 at->atom[ats[atHH]].m = at->atom[ats[atHH]].mB = 0;
1333 /* get dummy mass type */
1334 tpM = vsite_nm2type("MW", atype);
1335 /* make space for 1 mass: shift HH only */
1340 fprintf(stderr, "Inserting 1 dummy mass at %d\n", (*o2n)[i0] + 1);
1343 for (j = i0; j < at->nr; j++)
1345 (*o2n)[j] = j + *nadd;
1347 newx->resize(at->nr + *nadd);
1348 srenew(*newatom, at->nr + *nadd);
1349 srenew(*newatomname, at->nr + *nadd);
1350 srenew(*newvsite_type, at->nr + *nadd);
1351 srenew(*newcgnr, at->nr + *nadd);
1352 (*newatomname)[at->nr + *nadd - 1] = nullptr;
1354 /* Calc the dummy mass initial position */
1355 rvec_sub(x[ats[atHH]], x[ats[atOH]], r1);
1357 rvec_add(r1, x[ats[atHH]], (*newx)[atM]);
1359 strcpy(name, "MW1");
1360 (*newatomname)[atM] = put_symtab(symtab, name);
1361 (*newatom)[atM].m = (*newatom)[atM].mB = mM;
1362 (*newatom)[atM].q = (*newatom)[atM].qB = 0.0;
1363 (*newatom)[atM].type = (*newatom)[atM].typeB = tpM;
1364 (*newatom)[atM].ptype = eptAtom;
1365 (*newatom)[atM].resind = at->atom[i0].resind;
1366 (*newatom)[atM].elem[0] = 'M';
1367 (*newatom)[atM].elem[1] = '\0';
1368 (*newvsite_type)[atM] = NOTSET;
1369 (*newcgnr)[atM] = (*cgnr)[i0];
1370 /* renumber cgnr: */
1371 for (i = i0; i < at->nr; i++)
1376 (*vsite_type)[ats[atHH]] = F_VSITE2;
1378 /* assume we also want the COH angle constrained: */
1379 tmp1 = bond_cc * std::cos(0.5 * ANGLE_6RING) + dCGCE * std::sin(ANGLE_6RING * 0.5) + bond_co;
1380 dCGM = std::sqrt(cosrule(tmp1, vdist, angle_coh));
1381 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], add_shift + atM, dCGM);
1382 my_add_param(&(plist[F_CONSTRNC]), ats[atOH], add_shift + atM, vdist);
1384 add_vsite2_param(&plist[F_VSITE2], ats[atHH], ats[atOH], add_shift + atM, 1.0 / 2.0);
1388 static int gen_vsites_his(t_atoms* at,
1390 gmx::ArrayRef<InteractionsOfType> plist,
1393 gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
1396 real a, b, alpha, dCGCE1, dCGNE2;
1397 real sinalpha, cosalpha;
1398 real xcom, ycom, mtot;
1399 real mG, mrest, mCE1, mNE2;
1400 real b_CG_ND1, b_ND1_CE1, b_CE1_NE2, b_CG_CD2, b_CD2_NE2;
1401 real b_ND1_HD1, b_NE2_HE2, b_CE1_HE1, b_CD2_HD2;
1402 real a_CG_ND1_CE1, a_CG_CD2_NE2, a_ND1_CE1_NE2, a_CE1_NE2_CD2;
1403 real a_NE2_CE1_HE1, a_NE2_CD2_HD2, a_CE1_ND1_HD1, a_CE1_NE2_HE2;
1406 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1420 real x[atNR], y[atNR];
1422 /* CG, CE1 and NE2 stay, each gets part of the total mass,
1423 rest gets virtualized */
1424 /* check number of atoms, 3 hydrogens may be missing: */
1425 /* assert( nrfound >= atNR-3 || nrfound <= atNR );
1426 * Don't understand the above logic. Shouldn't it be && rather than || ???
1428 if ((nrfound < atNR - 3) || (nrfound > atNR))
1430 gmx_incons("Generating vsites for HIS");
1433 /* avoid warnings about uninitialized variables */
1434 b_ND1_HD1 = b_NE2_HE2 = b_CE1_HE1 = b_CD2_HD2 = a_NE2_CE1_HE1 = a_NE2_CD2_HD2 = a_CE1_ND1_HD1 =
1437 if (ats[atHD1] != NOTSET)
1439 if (ats[atHE2] != NOTSET)
1441 sprintf(resname, "HISH");
1445 sprintf(resname, "HISA");
1450 sprintf(resname, "HISB");
1453 /* Get geometry from database */
1454 b_CG_ND1 = get_ddb_bond(vsitetop, resname, "CG", "ND1");
1455 b_ND1_CE1 = get_ddb_bond(vsitetop, resname, "ND1", "CE1");
1456 b_CE1_NE2 = get_ddb_bond(vsitetop, resname, "CE1", "NE2");
1457 b_CG_CD2 = get_ddb_bond(vsitetop, resname, "CG", "CD2");
1458 b_CD2_NE2 = get_ddb_bond(vsitetop, resname, "CD2", "NE2");
1459 a_CG_ND1_CE1 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CG", "ND1", "CE1");
1460 a_CG_CD2_NE2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CG", "CD2", "NE2");
1461 a_ND1_CE1_NE2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "ND1", "CE1", "NE2");
1462 a_CE1_NE2_CD2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CE1", "NE2", "CD2");
1464 if (ats[atHD1] != NOTSET)
1466 b_ND1_HD1 = get_ddb_bond(vsitetop, resname, "ND1", "HD1");
1467 a_CE1_ND1_HD1 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CE1", "ND1", "HD1");
1469 if (ats[atHE2] != NOTSET)
1471 b_NE2_HE2 = get_ddb_bond(vsitetop, resname, "NE2", "HE2");
1472 a_CE1_NE2_HE2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CE1", "NE2", "HE2");
1474 if (ats[atHD2] != NOTSET)
1476 b_CD2_HD2 = get_ddb_bond(vsitetop, resname, "CD2", "HD2");
1477 a_NE2_CD2_HD2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "NE2", "CD2", "HD2");
1479 if (ats[atHE1] != NOTSET)
1481 b_CE1_HE1 = get_ddb_bond(vsitetop, resname, "CE1", "HE1");
1482 a_NE2_CE1_HE1 = DEG2RAD * get_ddb_angle(vsitetop, resname, "NE2", "CE1", "HE1");
1485 /* constraints between CG, CE1 and NE1 */
1486 dCGCE1 = std::sqrt(cosrule(b_CG_ND1, b_ND1_CE1, a_CG_ND1_CE1));
1487 dCGNE2 = std::sqrt(cosrule(b_CG_CD2, b_CD2_NE2, a_CG_CD2_NE2));
1489 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE1);
1490 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atNE2], dCGNE2);
1491 /* we already have a constraint CE1-NE2, so we don't add it again */
1493 /* calculate the positions in a local frame of reference.
1494 * The x-axis is the line from CG that makes a right angle
1495 * with the bond CE1-NE2, and the y-axis the bond CE1-NE2.
1497 /* First calculate the x-axis intersection with y-axis (=yCE1).
1498 * Get cos(angle CG-CE1-NE2) :
1500 cosalpha = acosrule(dCGNE2, dCGCE1, b_CE1_NE2);
1502 y[atCE1] = cosalpha * dCGCE1;
1504 y[atNE2] = y[atCE1] - b_CE1_NE2;
1505 sinalpha = std::sqrt(1 - cosalpha * cosalpha);
1506 x[atCG] = -sinalpha * dCGCE1;
1508 x[atHE1] = x[atHE2] = x[atHD1] = x[atHD2] = 0;
1509 y[atHE1] = y[atHE2] = y[atHD1] = y[atHD2] = 0;
1511 /* calculate ND1 and CD2 positions from CE1 and NE2 */
1513 x[atND1] = -b_ND1_CE1 * std::sin(a_ND1_CE1_NE2);
1514 y[atND1] = y[atCE1] - b_ND1_CE1 * std::cos(a_ND1_CE1_NE2);
1516 x[atCD2] = -b_CD2_NE2 * std::sin(a_CE1_NE2_CD2);
1517 y[atCD2] = y[atNE2] + b_CD2_NE2 * std::cos(a_CE1_NE2_CD2);
1519 /* And finally the hydrogen positions */
1520 if (ats[atHE1] != NOTSET)
1522 x[atHE1] = x[atCE1] + b_CE1_HE1 * std::sin(a_NE2_CE1_HE1);
1523 y[atHE1] = y[atCE1] - b_CE1_HE1 * std::cos(a_NE2_CE1_HE1);
1525 /* HD2 - first get (ccw) angle from (positive) y-axis */
1526 if (ats[atHD2] != NOTSET)
1528 alpha = a_CE1_NE2_CD2 + M_PI - a_NE2_CD2_HD2;
1529 x[atHD2] = x[atCD2] - b_CD2_HD2 * std::sin(alpha);
1530 y[atHD2] = y[atCD2] + b_CD2_HD2 * std::cos(alpha);
1532 if (ats[atHD1] != NOTSET)
1534 /* HD1 - first get (cw) angle from (positive) y-axis */
1535 alpha = a_ND1_CE1_NE2 + M_PI - a_CE1_ND1_HD1;
1536 x[atHD1] = x[atND1] - b_ND1_HD1 * std::sin(alpha);
1537 y[atHD1] = y[atND1] - b_ND1_HD1 * std::cos(alpha);
1539 if (ats[atHE2] != NOTSET)
1541 x[atHE2] = x[atNE2] + b_NE2_HE2 * std::sin(a_CE1_NE2_HE2);
1542 y[atHE2] = y[atNE2] + b_NE2_HE2 * std::cos(a_CE1_NE2_HE2);
1544 /* Have all coordinates now */
1546 /* calc center-of-mass; keep atoms CG, CE1, NE2 and
1547 * set the rest to vsite3
1549 mtot = xcom = ycom = 0;
1551 for (i = 0; i < atNR; i++)
1553 if (ats[i] != NOTSET)
1555 mtot += at->atom[ats[i]].m;
1556 xcom += x[i] * at->atom[ats[i]].m;
1557 ycom += y[i] * at->atom[ats[i]].m;
1558 if (i != atCG && i != atCE1 && i != atNE2)
1560 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
1561 (*vsite_type)[ats[i]] = F_VSITE3;
1566 if (nvsite + 3 != nrfound)
1568 gmx_incons("Generating vsites for HIS");
1574 /* distribute mass so that com stays the same */
1575 mG = xcom * mtot / x[atCG];
1577 mCE1 = (ycom - y[atNE2]) * mrest / (y[atCE1] - y[atNE2]);
1578 mNE2 = mrest - mCE1;
1580 at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = mG;
1581 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB = mCE1;
1582 at->atom[ats[atNE2]].m = at->atom[ats[atNE2]].mB = mNE2;
1585 if (ats[atHE1] != NOTSET)
1587 calc_vsite3_param(x[atHE1], y[atHE1], x[atCE1], y[atCE1], x[atNE2], y[atNE2], x[atCG],
1589 add_vsite3_param(&plist[F_VSITE3], ats[atHE1], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1592 if (ats[atHE2] != NOTSET)
1594 calc_vsite3_param(x[atHE2], y[atHE2], x[atNE2], y[atNE2], x[atCE1], y[atCE1], x[atCG],
1596 add_vsite3_param(&plist[F_VSITE3], ats[atHE2], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1600 calc_vsite3_param(x[atND1], y[atND1], x[atNE2], y[atNE2], x[atCE1], y[atCE1], x[atCG], y[atCG],
1602 add_vsite3_param(&plist[F_VSITE3], ats[atND1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1605 calc_vsite3_param(x[atCD2], y[atCD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2], x[atCG], y[atCG],
1607 add_vsite3_param(&plist[F_VSITE3], ats[atCD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1610 if (ats[atHD1] != NOTSET)
1612 calc_vsite3_param(x[atHD1], y[atHD1], x[atNE2], y[atNE2], x[atCE1], y[atCE1], x[atCG],
1614 add_vsite3_param(&plist[F_VSITE3], ats[atHD1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1617 if (ats[atHD2] != NOTSET)
1619 calc_vsite3_param(x[atHD2], y[atHD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2], x[atCG],
1621 add_vsite3_param(&plist[F_VSITE3], ats[atHD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1626 static bool is_vsite(int vsite_type)
1628 if (vsite_type == NOTSET)
1632 switch (abs(vsite_type))
1639 case F_VSITE4FDN: return TRUE;
1640 default: return FALSE;
1644 static char atomnamesuffix[] = "1234";
1646 void do_vsites(gmx::ArrayRef<const PreprocessResidue> rtpFFDB,
1647 PreprocessingAtomTypes* atype,
1650 std::vector<gmx::RVec>* x,
1651 gmx::ArrayRef<InteractionsOfType> plist,
1655 bool bVsiteAromatics,
1658 #define MAXATOMSPERRESIDUE 16
1659 int k, m, i0, ni0, whatres, add_shift, nvsite, nadd;
1661 int nrfound = 0, needed, nrbonds, nrHatoms, Heavy, nrheavies, tpM, tpHeavy;
1662 int Hatoms[4], heavies[4];
1663 bool bWARNING, bAddVsiteParam, bFirstWater;
1665 real mHtot, mtot, fact, fact2;
1666 rvec rpar, rperp, temp;
1667 char tpname[32], nexttpname[32];
1668 int * o2n, *newvsite_type, *newcgnr, ats[MAXATOMSPERRESIDUE];
1670 char*** newatomname;
1672 bool isN, planarN, bFound;
1674 /* if bVsiteAromatics=TRUE do_vsites will specifically convert atoms in
1675 PHE, TRP, TYR and HIS to a construction of virtual sites */
1684 const char* resnms[resNR] = { "PHE", "TRP", "TYR", "HIS" };
1685 /* Amber03 alternative names for termini */
1686 const char* resnmsN[resNR] = { "NPHE", "NTRP", "NTYR", "NHIS" };
1687 const char* resnmsC[resNR] = { "CPHE", "CTRP", "CTYR", "CHIS" };
1688 /* HIS can be known as HISH, HIS1, HISA, HID, HIE, HIP, etc. too */
1689 bool bPartial[resNR] = { FALSE, FALSE, FALSE, TRUE };
1690 /* the atnms for every residue MUST correspond to the enums in the
1691 gen_vsites_* (one for each residue) routines! */
1692 /* also the atom names in atnms MUST be in the same order as in the .rtp! */
1693 const char* atnms[resNR][MAXATOMSPERRESIDUE + 1] = {
1695 "CD1", "HD1", "CD2", "HD2", "CE1", "HE1", "CE2", "HE2", "CZ", "HZ", nullptr },
1697 "CG", "CD1", "HD1", "CD2", "NE1", "HE1", "CE2", "CE3", "HE3", "CZ2", "HZ2", "CZ3", "HZ3",
1698 "CH2", "HH2", nullptr },
1700 "CD1", "HD1", "CD2", "HD2", "CE1", "HE1", "CE2", "HE2", "CZ", "OH", "HH", nullptr },
1702 "ND1", "HD1", "CD2", "HD2", "CE1", "HE1", "NE2", "HE2", nullptr }
1707 printf("Searching for atoms to make virtual sites ...\n");
1708 fprintf(debug, "# # # VSITES # # #\n");
1711 std::vector<std::string> db = fflib_search_file_end(ffdir, ".vsd", FALSE);
1713 /* Container of CH3/NH3/NH2 configuration entries.
1714 * See comments in read_vsite_database. It isnt beautiful,
1715 * but it had to be fixed, and I dont even want to try to
1716 * maintain this part of the code...
1718 std::vector<VirtualSiteConfiguration> vsiteconflist;
1720 // TODO those have been deprecated and should be removed completely.
1721 /* Container of geometry (bond/angle) entries for
1722 * residues like PHE, TRP, TYR, HIS, etc., where we need
1723 * to know the geometry to construct vsite aromatics.
1724 * Note that equilibrium geometry isnt necessarily the same
1725 * as the individual bond and angle values given in the
1726 * force field (rings can be strained).
1728 std::vector<VirtualSiteTopology> vsitetop;
1729 for (const auto& filename : db)
1731 read_vsite_database(filename.c_str(), &vsiteconflist, &vsitetop);
1737 /* we need a marker for which atoms should *not* be renumbered afterwards */
1738 add_shift = 10 * at->nr;
1739 /* make arrays where masses can be inserted into */
1740 std::vector<gmx::RVec> newx(at->nr);
1741 snew(newatom, at->nr);
1742 snew(newatomname, at->nr);
1743 snew(newvsite_type, at->nr);
1744 snew(newcgnr, at->nr);
1745 /* make index array to tell where the atoms go to when masses are inserted */
1747 for (int i = 0; i < at->nr; i++)
1751 /* make index to tell which residues were already processed */
1752 std::vector<bool> bResProcessed(at->nres);
1756 /* generate vsite constructions */
1757 /* loop over all atoms */
1759 for (int i = 0; (i < at->nr); i++)
1761 if (at->atom[i].resind != resind)
1763 resind = at->atom[i].resind;
1765 const char* resnm = *(at->resinfo[resind].name);
1766 /* first check for aromatics to virtualize */
1767 /* don't waste our effort on DNA, water etc. */
1768 /* Only do the vsite aromatic stuff when we reach the
1769 * CA atom, since there might be an X2/X3 group on the
1770 * N-terminus that must be treated first.
1772 if (bVsiteAromatics && (strcmp(*(at->atomname[i]), "CA") == 0) && !bResProcessed[resind]
1773 && rt.namedResidueHasType(*(at->resinfo[resind].name), "Protein"))
1775 /* mark this residue */
1776 bResProcessed[resind] = TRUE;
1777 /* find out if this residue needs converting */
1779 for (int j = 0; j < resNR && whatres == NOTSET; j++)
1782 cmplength = bPartial[j] ? strlen(resnm) - 1 : strlen(resnm);
1784 bFound = ((gmx::equalCaseInsensitive(resnm, resnms[j], cmplength))
1785 || (gmx::equalCaseInsensitive(resnm, resnmsN[j], cmplength))
1786 || (gmx::equalCaseInsensitive(resnm, resnmsC[j], cmplength)));
1791 /* get atoms we will be needing for the conversion */
1793 for (k = 0; atnms[j][k]; k++)
1796 for (m = i; m < at->nr && at->atom[m].resind == resind && ats[k] == NOTSET; m++)
1798 if (gmx_strcasecmp(*(at->atomname[m]), atnms[j][k]) == 0)
1806 /* now k is number of atom names in atnms[j] */
1815 if (nrfound < needed)
1818 "not enough atoms found (%d, need %d) in "
1819 "residue %s %d while\n "
1820 "generating aromatics virtual site construction",
1821 nrfound, needed, resnm, at->resinfo[resind].nr);
1823 /* Advance overall atom counter */
1827 /* the enums for every residue MUST correspond to atnms[residue] */
1833 fprintf(stderr, "PHE at %d\n", o2n[ats[0]] + 1);
1835 nvsite += gen_vsites_phe(at, vsite_type, plist, nrfound, ats, vsitetop);
1840 fprintf(stderr, "TRP at %d\n", o2n[ats[0]] + 1);
1842 nvsite += gen_vsites_trp(atype, &newx, &newatom, &newatomname, &o2n,
1843 &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr, at,
1844 vsite_type, plist, nrfound, ats, add_shift, vsitetop);
1849 fprintf(stderr, "TYR at %d\n", o2n[ats[0]] + 1);
1851 nvsite += gen_vsites_tyr(atype, &newx, &newatom, &newatomname, &o2n,
1852 &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr, at,
1853 vsite_type, plist, nrfound, ats, add_shift, vsitetop);
1858 fprintf(stderr, "HIS at %d\n", o2n[ats[0]] + 1);
1860 nvsite += gen_vsites_his(at, vsite_type, plist, nrfound, ats, vsitetop);
1863 /* this means this residue won't be processed */
1865 default: gmx_fatal(FARGS, "DEATH HORROR in do_vsites (%s:%d)", __FILE__, __LINE__);
1866 } /* switch whatres */
1867 /* skip back to beginning of residue */
1868 while (i > 0 && at->atom[i - 1].resind == resind)
1872 } /* if bVsiteAromatics & is protein */
1874 /* now process the rest of the hydrogens */
1875 /* only process hydrogen atoms which are not already set */
1876 if (((*vsite_type)[i] == NOTSET) && is_hydrogen(*(at->atomname[i])))
1878 /* find heavy atom, count #bonds from it and #H atoms bound to it
1879 and return H atom numbers (Hatoms) and heavy atom numbers (heavies) */
1880 count_bonds(i, &plist[F_BONDS], at->atomname, &nrbonds, &nrHatoms, Hatoms, &Heavy,
1881 &nrheavies, heavies);
1882 /* get Heavy atom type */
1883 tpHeavy = get_atype(Heavy, at, rtpFFDB, &rt);
1884 strcpy(tpname, atype->atomNameFromAtomType(tpHeavy));
1887 bAddVsiteParam = TRUE;
1888 /* nested if's which check nrHatoms, nrbonds and atomname */
1893 case 2: /* -O-H */ (*vsite_type)[i] = F_BONDS; break;
1894 case 3: /* =CH-, -NH- or =NH+- */ (*vsite_type)[i] = F_VSITE3FD; break;
1895 case 4: /* --CH- (tert) */
1896 /* The old type 4FD had stability issues, so
1897 * all new constructs should use 4FDN
1899 (*vsite_type)[i] = F_VSITE4FDN;
1901 /* Check parity of heavy atoms from coordinates */
1906 rvec_sub((*x)[aj], (*x)[ai], tmpmat[0]);
1907 rvec_sub((*x)[ak], (*x)[ai], tmpmat[1]);
1908 rvec_sub((*x)[al], (*x)[ai], tmpmat[2]);
1910 if (det(tmpmat) > 0)
1918 default: /* nrbonds != 2, 3 or 4 */ bWARNING = TRUE;
1921 else if ((nrHatoms == 2) && (nrbonds == 2) && (at->atom[Heavy].atomnumber == 8))
1923 bAddVsiteParam = FALSE; /* this is water: skip these hydrogens */
1926 bFirstWater = FALSE;
1929 fprintf(debug, "Not converting hydrogens in water to virtual sites\n");
1933 else if ((nrHatoms == 2) && (nrbonds == 4))
1935 /* -CH2- , -NH2+- */
1936 (*vsite_type)[Hatoms[0]] = F_VSITE3OUT;
1937 (*vsite_type)[Hatoms[1]] = -F_VSITE3OUT;
1941 /* 2 or 3 hydrogen atom, with 3 or 4 bonds in total to the heavy atom.
1942 * If it is a nitrogen, first check if it is planar.
1944 isN = planarN = FALSE;
1945 if ((nrHatoms == 2) && ((*at->atomname[Heavy])[0] == 'N'))
1948 int j = nitrogen_is_planar(vsiteconflist, tpname);
1951 gmx_fatal(FARGS, "No vsite database NH2 entry for type %s\n", tpname);
1955 if ((nrHatoms == 2) && (nrbonds == 3) && (!isN || planarN))
1957 /* =CH2 or, if it is a nitrogen NH2, it is a planar one */
1958 (*vsite_type)[Hatoms[0]] = F_VSITE3FAD;
1959 (*vsite_type)[Hatoms[1]] = -F_VSITE3FAD;
1961 else if (((nrHatoms == 2) && (nrbonds == 3) && (isN && !planarN))
1962 || ((nrHatoms == 3) && (nrbonds == 4)))
1964 /* CH3, NH3 or non-planar NH2 group */
1965 int Hat_vsite_type[3] = { F_VSITE3, F_VSITE3OUT, F_VSITE3OUT };
1966 bool Hat_SwapParity[3] = { FALSE, TRUE, FALSE };
1970 fprintf(stderr, "-XH3 or nonplanar NH2 group at %d\n", i + 1);
1972 bAddVsiteParam = FALSE; /* we'll do this ourselves! */
1973 /* -NH2 (umbrella), -NH3+ or -CH3 */
1974 (*vsite_type)[Heavy] = F_VSITE3;
1975 for (int j = 0; j < nrHatoms; j++)
1977 (*vsite_type)[Hatoms[j]] = Hat_vsite_type[j];
1979 /* get dummy mass type from first char of heavy atom type (N or C) */
1982 atype->atomNameFromAtomType(get_atype(heavies[0], at, rtpFFDB, &rt)));
1983 std::string ch = get_dummymass_name(vsiteconflist, tpname, nexttpname);
1991 "Can't find dummy mass for type %s bonded to type %s in the "
1992 "virtual site database (.vsd files). Add it to the database!\n",
1993 tpname, nexttpname);
1998 "A dummy mass for type %s bonded to type %s is required, but "
1999 "no virtual site database (.vsd) files where found.\n",
2000 tpname, nexttpname);
2008 tpM = vsite_nm2type(name.c_str(), atype);
2009 /* make space for 2 masses: shift all atoms starting with 'Heavy' */
2015 fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, o2n[i0] + 1);
2018 for (int j = i0; j < at->nr; j++)
2023 newx.resize(at->nr + nadd);
2024 srenew(newatom, at->nr + nadd);
2025 srenew(newatomname, at->nr + nadd);
2026 srenew(newvsite_type, at->nr + nadd);
2027 srenew(newcgnr, at->nr + nadd);
2029 for (int j = 0; j < NMASS; j++)
2031 newatomname[at->nr + nadd - 1 - j] = nullptr;
2034 /* calculate starting position for the masses */
2036 /* get atom masses, and set Heavy and Hatoms mass to zero */
2037 for (int j = 0; j < nrHatoms; j++)
2039 mHtot += get_amass(Hatoms[j], at, rtpFFDB, &rt);
2040 at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
2042 mtot = mHtot + get_amass(Heavy, at, rtpFFDB, &rt);
2043 at->atom[Heavy].m = at->atom[Heavy].mB = 0;
2048 fact2 = mHtot / mtot;
2049 fact = std::sqrt(fact2);
2050 /* generate vectors parallel and perpendicular to rotational axis:
2051 * rpar = Heavy -> Hcom
2052 * rperp = Hcom -> H1 */
2054 for (int j = 0; j < nrHatoms; j++)
2056 rvec_inc(rpar, (*x)[Hatoms[j]]);
2058 svmul(1.0 / nrHatoms, rpar, rpar); /* rpar = ( H1+H2+H3 ) / 3 */
2059 rvec_dec(rpar, (*x)[Heavy]); /* - Heavy */
2060 rvec_sub((*x)[Hatoms[0]], (*x)[Heavy], rperp);
2061 rvec_dec(rperp, rpar); /* rperp = H1 - Heavy - rpar */
2062 /* calc mass positions */
2063 svmul(fact2, rpar, temp);
2064 for (int j = 0; (j < NMASS); j++) /* xM = xN + fact2 * rpar +/- fact * rperp */
2066 rvec_add((*x)[Heavy], temp, newx[ni0 + j]);
2068 svmul(fact, rperp, temp);
2069 rvec_inc(newx[ni0], temp);
2070 rvec_dec(newx[ni0 + 1], temp);
2071 /* set atom parameters for the masses */
2072 for (int j = 0; (j < NMASS); j++)
2074 /* make name: "M??#" or "M?#" (? is atomname, # is number) */
2077 for (k = 0; (*at->atomname[Heavy])[k] && (k < NMASS); k++)
2079 name[k + 1] = (*at->atomname[Heavy])[k];
2081 name[k + 1] = atomnamesuffix[j];
2083 newatomname[ni0 + j] = put_symtab(symtab, name.c_str());
2084 newatom[ni0 + j].m = newatom[ni0 + j].mB = mtot / NMASS;
2085 newatom[ni0 + j].q = newatom[ni0 + j].qB = 0.0;
2086 newatom[ni0 + j].type = newatom[ni0 + j].typeB = tpM;
2087 newatom[ni0 + j].ptype = eptAtom;
2088 newatom[ni0 + j].resind = at->atom[i0].resind;
2089 newatom[ni0 + j].elem[0] = 'M';
2090 newatom[ni0 + j].elem[1] = '\0';
2091 newvsite_type[ni0 + j] = NOTSET;
2092 newcgnr[ni0 + j] = (*cgnr)[i0];
2094 /* add constraints between dummy masses and to heavies[0] */
2095 /* 'add_shift' says which atoms won't be renumbered afterwards */
2096 my_add_param(&(plist[F_CONSTRNC]), heavies[0], add_shift + ni0, NOTSET);
2097 my_add_param(&(plist[F_CONSTRNC]), heavies[0], add_shift + ni0 + 1, NOTSET);
2098 my_add_param(&(plist[F_CONSTRNC]), add_shift + ni0, add_shift + ni0 + 1, NOTSET);
2100 /* generate Heavy, H1, H2 and H3 from M1, M2 and heavies[0] */
2101 /* note that vsite_type cannot be NOTSET, because we just set it */
2102 add_vsite3_atoms(&plist[(*vsite_type)[Heavy]], Heavy, heavies[0],
2103 add_shift + ni0, add_shift + ni0 + 1, FALSE);
2104 for (int j = 0; j < nrHatoms; j++)
2106 add_vsite3_atoms(&plist[(*vsite_type)[Hatoms[j]]], Hatoms[j], heavies[0],
2107 add_shift + ni0, add_shift + ni0 + 1, Hat_SwapParity[j]);
2119 "Cannot convert atom %d %s (bound to a heavy atom "
2121 " %d bonds and %d bound hydrogens atoms) to virtual site\n",
2122 i + 1, *(at->atomname[i]), tpname, nrbonds, nrHatoms);
2126 /* add vsite parameters to topology,
2127 also get rid of negative vsite_types */
2128 add_vsites(plist, (*vsite_type), Heavy, nrHatoms, Hatoms, nrheavies, heavies);
2129 /* transfer mass of virtual site to Heavy atom */
2130 for (int j = 0; j < nrHatoms; j++)
2132 if (is_vsite((*vsite_type)[Hatoms[j]]))
2134 at->atom[Heavy].m += at->atom[Hatoms[j]].m;
2135 at->atom[Heavy].mB = at->atom[Heavy].m;
2136 at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
2143 fprintf(debug, "atom %d: ", o2n[i] + 1);
2144 print_bonds(debug, o2n, nrHatoms, Hatoms, Heavy, nrheavies, heavies);
2146 } /* if vsite NOTSET & is hydrogen */
2148 } /* for i < at->nr */
2152 fprintf(debug, "Before inserting new atoms:\n");
2153 for (int i = 0; i < at->nr; i++)
2155 fprintf(debug, "%4d %4d %4s %4d %4s %6d %-10s\n", i + 1, o2n[i] + 1,
2156 at->atomname[i] ? *(at->atomname[i]) : "(NULL)", at->resinfo[at->atom[i].resind].nr,
2157 at->resinfo[at->atom[i].resind].name ? *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
2159 ((*vsite_type)[i] == NOTSET) ? "NOTSET" : interaction_function[(*vsite_type)[i]].name);
2161 fprintf(debug, "new atoms to be inserted:\n");
2162 for (int i = 0; i < at->nr + nadd; i++)
2166 fprintf(debug, "%4d %4s %4d %6d %-10s\n", i + 1,
2167 newatomname[i] ? *(newatomname[i]) : "(NULL)", newatom[i].resind, newcgnr[i],
2168 (newvsite_type[i] == NOTSET) ? "NOTSET"
2169 : interaction_function[newvsite_type[i]].name);
2174 /* add all original atoms to the new arrays, using o2n index array */
2175 for (int i = 0; i < at->nr; i++)
2177 newatomname[o2n[i]] = at->atomname[i];
2178 newatom[o2n[i]] = at->atom[i];
2179 newvsite_type[o2n[i]] = (*vsite_type)[i];
2180 newcgnr[o2n[i]] = (*cgnr)[i];
2181 copy_rvec((*x)[i], newx[o2n[i]]);
2183 /* throw away old atoms */
2185 sfree(at->atomname);
2188 /* put in the new ones */
2191 at->atomname = newatomname;
2192 *vsite_type = newvsite_type;
2195 if (at->nr > add_shift)
2198 "Added impossible amount of dummy masses "
2199 "(%d on a total of %d atoms)\n",
2200 nadd, at->nr - nadd);
2205 fprintf(debug, "After inserting new atoms:\n");
2206 for (int i = 0; i < at->nr; i++)
2208 fprintf(debug, "%4d %4s %4d %4s %6d %-10s\n", i + 1,
2209 at->atomname[i] ? *(at->atomname[i]) : "(NULL)", at->resinfo[at->atom[i].resind].nr,
2210 at->resinfo[at->atom[i].resind].name ? *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
2212 ((*vsite_type)[i] == NOTSET) ? "NOTSET" : interaction_function[(*vsite_type)[i]].name);
2216 /* now renumber all the interactions because of the added atoms */
2217 for (int ftype = 0; ftype < F_NRE; ftype++)
2219 InteractionsOfType* params = &(plist[ftype]);
2222 fprintf(debug, "Renumbering %zu %s\n", params->size(), interaction_function[ftype].longname);
2224 /* Horrible hacks needed here to get this to work */
2225 for (auto parm = params->interactionTypes.begin(); parm != params->interactionTypes.end(); parm++)
2227 gmx::ArrayRef<const int> atomNumbers(parm->atoms());
2228 std::vector<int> newAtomNumber;
2229 for (int j = 0; j < NRAL(ftype); j++)
2231 if (atomNumbers[j] >= add_shift)
2235 fprintf(debug, " [%d -> %d]", atomNumbers[j], atomNumbers[j] - add_shift);
2237 newAtomNumber.emplace_back(atomNumbers[j] - add_shift);
2243 fprintf(debug, " [%d -> %d]", atomNumbers[j], o2n[atomNumbers[j]]);
2245 newAtomNumber.emplace_back(o2n[atomNumbers[j]]);
2248 *parm = InteractionOfType(newAtomNumber, parm->forceParam(), parm->interactionTypeName());
2251 fprintf(debug, "\n");
2255 /* sort constraint parameters */
2256 InteractionsOfType* params = &(plist[F_CONSTRNC]);
2257 for (auto& type : params->interactionTypes)
2265 /* tell the user what we did */
2266 fprintf(stderr, "Marked %d virtual sites\n", nvsite);
2267 fprintf(stderr, "Added %d dummy masses\n", nadd);
2268 fprintf(stderr, "Added %zu new constraints\n", plist[F_CONSTRNC].size());
2271 void do_h_mass(InteractionsOfType* psb, int vsite_type[], t_atoms* at, real mHmult, bool bDeuterate)
2273 /* loop over all atoms */
2274 for (int i = 0; i < at->nr; i++)
2276 /* adjust masses if i is hydrogen and not a virtual site */
2277 if (!is_vsite(vsite_type[i]) && is_hydrogen(*(at->atomname[i])))
2279 /* find bonded heavy atom */
2281 for (auto parm = psb->interactionTypes.begin();
2282 (parm != psb->interactionTypes.end()) && (a == NOTSET); parm++)
2284 /* if other atom is not a virtual site, it is the one we want */
2285 if ((parm->ai() == i) && !is_vsite(vsite_type[parm->aj()]))
2289 else if ((parm->aj() == i) && !is_vsite(vsite_type[parm->ai()]))
2296 gmx_fatal(FARGS, "Unbound hydrogen atom (%d) found while adjusting mass", i + 1);
2299 /* adjust mass of i (hydrogen) with mHmult
2300 and correct mass of a (bonded atom) with same amount */
2303 at->atom[a].m -= (mHmult - 1.0) * at->atom[i].m;
2304 at->atom[a].mB -= (mHmult - 1.0) * at->atom[i].m;
2306 at->atom[i].m *= mHmult;
2307 at->atom[i].mB *= mHmult;