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,2018,2019, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 #include "gen_vsite.h"
51 #include "gromacs/fileio/pdbio.h"
52 #include "gromacs/gmxpreprocess/add_par.h"
53 #include "gromacs/gmxpreprocess/fflibutil.h"
54 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
55 #include "gromacs/gmxpreprocess/grompp_impl.h"
56 #include "gromacs/gmxpreprocess/notset.h"
57 #include "gromacs/gmxpreprocess/toputil.h"
58 #include "gromacs/math/functions.h"
59 #include "gromacs/math/units.h"
60 #include "gromacs/math/vec.h"
61 #include "gromacs/mdtypes/md_enums.h"
62 #include "gromacs/topology/ifunc.h"
63 #include "gromacs/topology/residuetypes.h"
64 #include "gromacs/topology/symtab.h"
65 #include "gromacs/utility/basedefinitions.h"
66 #include "gromacs/utility/cstringutil.h"
67 #include "gromacs/utility/fatalerror.h"
68 #include "gromacs/utility/futil.h"
69 #include "gromacs/utility/real.h"
70 #include "gromacs/utility/smalloc.h"
72 #include "hackblock.h"
76 #define OPENDIR '[' /* starting sign for directive */
77 #define CLOSEDIR ']' /* ending sign for directive */
79 /*! \libinternal \brief
80 * The configuration describing a virtual site.
82 struct VirtualSiteConfiguration
85 * Explicit constructor.
87 * \param[in] type Atomtype for vsite configuration.
88 * \param[in] planar Is the input conf planar.
89 * \param[in] nhyd How many hydrogens are in the configuration.
90 * \param[in] nextheavy Type of bonded heavy atom.
91 * \param[in] dummy What kind of dummy is used in the vsite.
93 explicit VirtualSiteConfiguration(const std::string &type, bool planar,
94 int nhyd, const std::string &nextheavy, const std::string &dummy)
95 : atomtype(type), isplanar(planar), nHydrogens(nhyd), nextHeavyType(nextheavy),
98 //! Type for the XH3/XH2 atom.
100 /*! \brief Is the configuration planar?
102 * If true, the atomtype above and the three connected
103 * ones are in a planar geometry. The two next entries
104 * are undefined in that case.
106 bool isplanar = false;
107 //!cnumber of connected hydrogens.
109 //! Type for the heavy atom bonded to XH2/XH3.
110 std::string nextHeavyType;
111 //! The type of MNH* or MCH3* dummy mass to use.
112 std::string dummyMass;
116 /*!\libinternal \brief
117 * Virtual site topology datastructure.
119 * Structure to represent average bond and angles values in vsite aromatic
120 * residues. Note that these are NOT necessarily the bonds and angles from the
121 * forcefield; many forcefields (like Amber, OPLS) have some inherent strain in
122 * 5-rings (i.e. the sum of angles is !=540, but impropers keep it planar)
124 struct VirtualSiteTopology
127 * Explicit constructor
129 * \param[in] name Residue name.
131 explicit VirtualSiteTopology(const std::string &name) : resname(name)
135 //! Helper struct for single bond in virtual site.
136 struct VirtualSiteBond
139 * Explicit constructor
141 * \param[in] a1 First atom name.
142 * \param[in] a2 Second atom name.
143 * \param[in] v Value for distance.
145 VirtualSiteBond(const std::string &a1, const std::string &a2, real v) :
146 atom1(a1), atom2(a2), value(v)
152 //! Distance value between atoms.
155 //! Container of all bonds in virtual site.
156 std::vector<VirtualSiteBond> bond;
157 //! Helper struct for single angle in virtual site.
158 struct VirtualSiteAngle
161 * Explicit constructor
163 * \param[in] a1 First atom name.
164 * \param[in] a2 Second atom name.
165 * \param[in] a3 Third atom name.
166 * \param[in] v Value for angle.
168 VirtualSiteAngle(const std::string &a1, const std::string &a2, const std::string &a3, real v) :
169 atom1(a1), atom2(a2), atom3(a3), value(v)
180 //! Container for all angles in virtual site.
181 std::vector<VirtualSiteAngle> angle;
186 DDB_CH3, DDB_NH3, DDB_NH2, DDB_PHE, DDB_TYR,
187 DDB_TRP, DDB_HISA, DDB_HISB, DDB_HISH, DDB_DIR_NR
190 typedef char t_dirname[STRLEN];
192 static const t_dirname ddb_dirnames[DDB_DIR_NR] = {
204 static int ddb_name2dir(char *name)
206 /* Translate a directive name to the number of the directive.
207 * HID/HIE/HIP names are translated to the ones we use in Gromacs.
214 for (i = 0; i < DDB_DIR_NR && index < 0; i++)
216 if (!gmx_strcasecmp(name, ddb_dirnames[i]))
226 static void read_vsite_database(const char *ddbname,
227 std::vector<VirtualSiteConfiguration> *vsiteconflist,
228 std::vector<VirtualSiteTopology> *vsitetoplist)
230 /* This routine is a quick hack to fix the problem with hardcoded atomtypes
231 * and aromatic vsite parameters by reading them from a ff???.vsd file.
233 * The file can contain sections [ NH3 ], [ CH3 ], [ NH2 ], and ring residue names.
234 * For the NH3 and CH3 section each line has three fields. The first is the atomtype
235 * (nb: not bonded type) of the N/C atom to be replaced, the second field is
236 * the type of the next heavy atom it is bonded to, and the third field the type
237 * of dummy mass that will be used for this group.
239 * If the NH2 group planar (sp2 N) a different vsite construct is used, so in this
240 * case the second field should just be the word planar.
247 char s1[MAXNAME], s2[MAXNAME], s3[MAXNAME], s4[MAXNAME];
249 gmx::FilePtr ddb = gmx::openLibraryFile(ddbname);
253 while (fgets2(pline, STRLEN-2, ddb.get()) != nullptr)
255 strip_comment(pline);
257 if (strlen(pline) > 0)
259 if (pline[0] == OPENDIR)
261 strncpy(dirstr, pline+1, STRLEN-2);
262 if ((ch = strchr (dirstr, CLOSEDIR)) != nullptr)
268 if (!gmx_strcasecmp(dirstr, "HID") ||
269 !gmx_strcasecmp(dirstr, "HISD"))
271 sprintf(dirstr, "HISA");
273 else if (!gmx_strcasecmp(dirstr, "HIE") ||
274 !gmx_strcasecmp(dirstr, "HISE"))
276 sprintf(dirstr, "HISB");
278 else if (!gmx_strcasecmp(dirstr, "HIP"))
280 sprintf(dirstr, "HISH");
283 curdir = ddb_name2dir(dirstr);
286 gmx_fatal(FARGS, "Invalid directive %s in vsite database %s",
295 gmx_fatal(FARGS, "First entry in vsite database must be a directive.\n");
300 int numberOfSites = sscanf(pline, "%s%s%s", s1, s2, s3);
301 std::string s1String = s1;
302 std::string s2String = s2;
303 std::string s3String = s3;
304 if (numberOfSites < 3 && gmx::equalCaseInsensitive(s2String, "planar"))
306 VirtualSiteConfiguration newVsiteConf(s1String, true, 2, "0", "0");
307 vsiteconflist->push_back(newVsiteConf);
309 else if (numberOfSites == 3)
311 VirtualSiteConfiguration newVsiteConf(s1String, false, -1, s2String, s3String);
312 if (curdir == DDB_NH2)
314 newVsiteConf.nHydrogens = 2;
318 newVsiteConf.nHydrogens = 3;
320 vsiteconflist->push_back(newVsiteConf);
324 gmx_fatal(FARGS, "Not enough directives in vsite database line: %s\n", pline);
335 const auto found = std::find_if(vsitetoplist->begin(), vsitetoplist->end(),
336 [&dirstr](const auto &entry)
337 { return gmx::equalCaseInsensitive(dirstr, entry.resname); });
338 /* Allocate a new topology entry if this is a new residue */
339 if (found == vsitetoplist->end())
341 vsitetoplist->push_back(VirtualSiteTopology(dirstr));
343 int numberOfSites = sscanf(pline, "%s%s%s%s", s1, s2, s3, s4);
344 std::string s1String = s1;
345 std::string s2String = s2;
346 std::string s3String = s3;
348 if (numberOfSites == 3)
351 vsitetoplist->back().bond.emplace_back(s1String, s2String, strtod(s3, nullptr));
353 else if (numberOfSites == 4)
356 vsitetoplist->back().angle.emplace_back(s1String, s2String, s3String, strtod(s4, nullptr));
361 gmx_fatal(FARGS, "Need 3 or 4 values to specify bond/angle values in %s: %s\n", ddbname, pline);
366 gmx_fatal(FARGS, "Didnt find a case for directive %s in read_vsite_database\n", dirstr);
373 static int nitrogen_is_planar(gmx::ArrayRef<const VirtualSiteConfiguration> vsiteconflist,
374 const std::string &atomtype)
376 /* Return 1 if atomtype exists in database list and is planar, 0 if not,
377 * and -1 if not found.
380 const auto found = std::find_if(vsiteconflist.begin(), vsiteconflist.end(),
381 [&atomtype](const auto &entry)
382 { return (gmx::equalCaseInsensitive(entry.atomtype, atomtype) &&
383 entry.nHydrogens == 2); });
384 if (found != vsiteconflist.end())
386 res = static_cast<int>(found->isplanar);
396 static const std::string get_dummymass_name(gmx::ArrayRef<const VirtualSiteConfiguration> vsiteconflist,
397 const std::string &atom, const std::string &nextheavy)
399 /* Return the dummy mass name if found, or NULL if not set in ddb database */
400 const auto found = std::find_if(vsiteconflist.begin(), vsiteconflist.end(),
401 [&atom, &nextheavy](const auto &entry)
402 { return (gmx::equalCaseInsensitive(atom, entry.atomtype) &&
403 gmx::equalCaseInsensitive(nextheavy, entry.nextHeavyType)); });
404 if (found != vsiteconflist.end())
406 return found->dummyMass;
416 static real get_ddb_bond(gmx::ArrayRef<const VirtualSiteTopology> vsitetop,
417 const std::string &res,
418 const std::string &atom1,
419 const std::string &atom2)
421 const auto found = std::find_if(vsitetop.begin(), vsitetop.end(),
422 [&res](const auto &entry)
423 { return gmx::equalCaseInsensitive(res, entry.resname); });
425 if (found == vsitetop.end())
427 gmx_fatal(FARGS, "No vsite information for residue %s found in vsite database.\n", res.c_str());
429 const auto foundBond = std::find_if(found->bond.begin(), found->bond.end(),
430 [&atom1, &atom2](const auto &entry)
431 { return ((atom1 == entry.atom1 && atom2 == entry.atom2) ||
432 (atom1 == entry.atom2 && atom2 == entry.atom1)); });
433 if (foundBond == found->bond.end())
435 gmx_fatal(FARGS, "Couldnt find bond %s-%s for residue %s in vsite database.\n", atom1.c_str(), atom2.c_str(), res.c_str());
438 return foundBond->value;
442 static real get_ddb_angle(gmx::ArrayRef<const VirtualSiteTopology> vsitetop,
443 const std::string &res,
444 const std::string &atom1,
445 const std::string &atom2,
446 const std::string &atom3)
448 const auto found = std::find_if(vsitetop.begin(), vsitetop.end(),
449 [&res](const auto &entry)
450 { return gmx::equalCaseInsensitive(res, entry.resname); });
452 if (found == vsitetop.end())
454 gmx_fatal(FARGS, "No vsite information for residue %s found in vsite database.\n", res.c_str());
456 const auto foundAngle = std::find_if(found->angle.begin(), found->angle.end(),
457 [&atom1, &atom2, &atom3](const auto &entry)
458 { return ((atom1 == entry.atom1 && atom2 == entry.atom2 && atom3 == entry.atom3) ||
459 (atom1 == entry.atom3 && atom2 == entry.atom2 && atom3 == entry.atom1) ||
460 (atom1 == entry.atom2 && atom2 == entry.atom1 && atom3 == entry.atom3) ||
461 (atom1 == entry.atom3 && atom2 == entry.atom1 && atom3 == entry.atom2)); });
463 if (foundAngle == found->angle.end())
465 gmx_fatal(FARGS, "Couldnt find angle %s-%s-%s for residue %s in vsite database.\n", atom1.c_str(), atom2.c_str(), atom3.c_str(), res.c_str());
468 return foundAngle->value;
472 static void count_bonds(int atom, InteractionTypeParameters *psb, char ***atomname,
473 int *nrbonds, int *nrHatoms, int Hatoms[], int *Heavy,
474 int *nrheavies, int heavies[])
476 int heavy, other, nrb, nrH, nrhv;
478 /* find heavy atom bound to this hydrogen */
480 for (int i = 0; (i < psb->nr) && (heavy == NOTSET); i++)
482 if (psb->param[i].ai() == atom)
484 heavy = psb->param[i].aj();
486 else if (psb->param[i].aj() == atom)
488 heavy = psb->param[i].ai();
493 gmx_fatal(FARGS, "unbound hydrogen atom %d", atom+1);
495 /* find all atoms bound to heavy atom */
500 for (int i = 0; i < psb->nr; i++)
502 if (psb->param[i].ai() == heavy)
504 other = psb->param[i].aj();
506 else if (psb->param[i].aj() == heavy)
508 other = psb->param[i].ai();
513 if (is_hydrogen(*(atomname[other])))
520 heavies[nrhv] = other;
532 static void print_bonds(FILE *fp, int o2n[],
533 int nrHatoms, const int Hatoms[], int Heavy,
534 int nrheavies, const int heavies[])
538 fprintf(fp, "Found: %d Hatoms: ", nrHatoms);
539 for (i = 0; i < nrHatoms; i++)
541 fprintf(fp, " %d", o2n[Hatoms[i]]+1);
543 fprintf(fp, "; %d Heavy atoms: %d", nrheavies+1, o2n[Heavy]+1);
544 for (i = 0; i < nrheavies; i++)
546 fprintf(fp, " %d", o2n[heavies[i]]+1);
551 static int get_atype(int atom, t_atoms *at, gmx::ArrayRef<const PreprocessResidue> rtpFFDB,
557 if (at->atom[atom].m != 0.0f)
559 type = at->atom[atom].type;
563 /* get type from rtpFFDB */
564 auto localPpResidue = getDatabaseEntry(*(at->resinfo[at->atom[atom].resind].name), rtpFFDB);
565 bNterm = rt->namedResidueHasType(*(at->resinfo[at->atom[atom].resind].name), "Protein") &&
566 (at->atom[atom].resind == 0);
567 int j = search_jtype(*localPpResidue, *(at->atomname[atom]), bNterm);
568 type = localPpResidue->atom[j].type;
573 static int vsite_nm2type(const char *name, PreprocessingAtomTypes *atype)
577 tp = atype->atomTypeFromName(name);
580 gmx_fatal(FARGS, "Dummy mass type (%s) not found in atom type database",
587 static real get_amass(int atom, t_atoms *at, gmx::ArrayRef<const PreprocessResidue> rtpFFDB,
593 if (at->atom[atom].m != 0.0f)
595 mass = at->atom[atom].m;
599 /* get mass from rtpFFDB */
600 auto localPpResidue = getDatabaseEntry(*(at->resinfo[at->atom[atom].resind].name), rtpFFDB);
601 bNterm = rt->namedResidueHasType(*(at->resinfo[at->atom[atom].resind].name), "Protein") &&
602 (at->atom[atom].resind == 0);
603 int j = search_jtype(*localPpResidue, *(at->atomname[atom]), bNterm);
604 mass = localPpResidue->atom[j].m;
609 static void my_add_param(InteractionTypeParameters *plist, int ai, int aj, real b)
611 static real c[MAXFORCEPARAM] =
612 { NOTSET, NOTSET, NOTSET, NOTSET, NOTSET, NOTSET };
615 add_param(plist, ai, aj, c, nullptr);
618 static void add_vsites(gmx::ArrayRef<InteractionTypeParameters> plist, int vsite_type[],
619 int Heavy, int nrHatoms, int Hatoms[],
620 int nrheavies, int heavies[])
622 int other, moreheavy;
624 for (int i = 0; i < nrHatoms; i++)
626 int ftype = vsite_type[Hatoms[i]];
627 /* Errors in setting the vsite_type should really be caugth earlier,
628 * because here it's not possible to print any useful error message.
629 * But it's still better to print a message than to segfault.
633 gmx_incons("Undetected error in setting up virtual sites");
635 bool bSwapParity = (ftype < 0);
636 vsite_type[Hatoms[i]] = ftype = abs(ftype);
637 if (ftype == F_BONDS)
639 if ( (nrheavies != 1) && (nrHatoms != 1) )
641 gmx_fatal(FARGS, "cannot make constraint in add_vsites for %d heavy "
642 "atoms and %d hydrogen atoms", nrheavies, nrHatoms);
644 my_add_param(&(plist[F_CONSTRNC]), Hatoms[i], heavies[0], NOTSET);
655 gmx_fatal(FARGS, "Not enough heavy atoms (%d) for %s (min 3)",
657 interaction_function[vsite_type[Hatoms[i]]].name);
659 add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], heavies[1],
666 moreheavy = heavies[1];
670 /* find more heavy atoms */
671 other = moreheavy = NOTSET;
672 for (int j = 0; (j < plist[F_BONDS].nr) && (moreheavy == NOTSET); j++)
674 if (plist[F_BONDS].param[j].ai() == heavies[0])
676 other = plist[F_BONDS].param[j].aj();
678 else if (plist[F_BONDS].param[j].aj() == heavies[0])
680 other = plist[F_BONDS].param[j].ai();
682 if ( (other != NOTSET) && (other != Heavy) )
687 if (moreheavy == NOTSET)
689 gmx_fatal(FARGS, "Unbound molecule part %d-%d", Heavy+1, Hatoms[0]+1);
692 add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], moreheavy,
700 gmx_fatal(FARGS, "Not enough heavy atoms (%d) for %s (min 4)",
702 interaction_function[vsite_type[Hatoms[i]]].name);
704 add_vsite4_atoms(&plist[ftype],
705 Hatoms[0], Heavy, heavies[0], heavies[1], heavies[2]);
709 gmx_fatal(FARGS, "can't use add_vsites for interaction function %s",
710 interaction_function[vsite_type[Hatoms[i]]].name);
716 #define ANGLE_6RING (DEG2RAD*120)
718 /* cosine rule: a^2 = b^2 + c^2 - 2 b c cos(alpha) */
719 /* get a^2 when a, b and alpha are given: */
720 #define cosrule(b, c, alpha) ( gmx::square(b) + gmx::square(c) - 2*(b)*(c)*std::cos(alpha) )
721 /* get cos(alpha) when a, b and c are given: */
722 #define acosrule(a, b, c) ( (gmx::square(b)+gmx::square(c)-gmx::square(a))/(2*(b)*(c)) )
724 static int gen_vsites_6ring(t_atoms *at, int *vsite_type[], gmx::ArrayRef<InteractionTypeParameters> plist,
725 int nrfound, int *ats, real bond_cc, real bond_ch,
726 real xcom, bool bDoZ)
728 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
730 atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
735 real a, b, dCGCE, tmp1, tmp2, mtot, mG, mrest;
737 /* CG, CE1 and CE2 stay and each get a part of the total mass,
738 * so the c-o-m stays the same.
745 gmx_incons("Generating vsites on 6-rings");
749 /* constraints between CG, CE1 and CE2: */
750 dCGCE = std::sqrt( cosrule(bond_cc, bond_cc, ANGLE_6RING) );
751 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE);
752 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE2], dCGCE);
753 my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atCE2], dCGCE);
755 /* rest will be vsite3 */
758 for (i = 0; i < (bDoZ ? atNR : atHZ); i++)
760 mtot += at->atom[ats[i]].m;
761 if (i != atCG && i != atCE1 && i != atCE2 && (bDoZ || (i != atHZ && i != atCZ) ) )
763 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
764 (*vsite_type)[ats[i]] = F_VSITE3;
768 /* Distribute mass so center-of-mass stays the same.
769 * The center-of-mass in the call is defined with x=0 at
770 * the CE1-CE2 bond and y=0 at the line from CG to the middle of CE1-CE2 bond.
772 xCG = -bond_cc+bond_cc*std::cos(ANGLE_6RING);
774 mG = at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = xcom*mtot/xCG;
776 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB =
777 at->atom[ats[atCE2]].m = at->atom[ats[atCE2]].mB = mrest / 2;
779 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
780 tmp1 = dCGCE*std::sin(ANGLE_6RING*0.5);
781 tmp2 = bond_cc*std::cos(0.5*ANGLE_6RING) + tmp1;
783 a = b = -bond_ch / tmp1;
785 add_vsite3_param(&plist[F_VSITE3],
786 ats[atHE1], ats[atCE1], ats[atCE2], ats[atCG], a, b);
787 add_vsite3_param(&plist[F_VSITE3],
788 ats[atHE2], ats[atCE2], ats[atCE1], ats[atCG], a, b);
789 /* CD1, CD2 and CZ: */
791 add_vsite3_param(&plist[F_VSITE3],
792 ats[atCD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
793 add_vsite3_param(&plist[F_VSITE3],
794 ats[atCD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
797 add_vsite3_param(&plist[F_VSITE3],
798 ats[atCZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
800 /* HD1, HD2 and HZ: */
801 a = b = ( bond_ch + tmp2 ) / tmp1;
802 add_vsite3_param(&plist[F_VSITE3],
803 ats[atHD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
804 add_vsite3_param(&plist[F_VSITE3],
805 ats[atHD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
808 add_vsite3_param(&plist[F_VSITE3],
809 ats[atHZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
815 static int gen_vsites_phe(t_atoms *at, int *vsite_type[], gmx::ArrayRef<InteractionTypeParameters> plist,
816 int nrfound, int *ats, gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
818 real bond_cc, bond_ch;
821 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
823 atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
827 /* Aromatic rings have 6-fold symmetry, so we only need one bond length.
828 * (angle is always 120 degrees).
830 bond_cc = get_ddb_bond(vsitetop, "PHE", "CD1", "CE1");
831 bond_ch = get_ddb_bond(vsitetop, "PHE", "CD1", "HD1");
833 x[atCG] = -bond_cc+bond_cc*std::cos(ANGLE_6RING);
835 x[atHD1] = x[atCD1]+bond_ch*std::cos(ANGLE_6RING);
837 x[atHE1] = x[atCE1]-bond_ch*std::cos(ANGLE_6RING);
842 x[atCZ] = bond_cc*std::cos(0.5*ANGLE_6RING);
843 x[atHZ] = x[atCZ]+bond_ch;
846 for (i = 0; i < atNR; i++)
848 xcom += x[i]*at->atom[ats[i]].m;
849 mtot += at->atom[ats[i]].m;
853 return gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, TRUE);
856 static void calc_vsite3_param(real xd, real yd, real xi, real yi, real xj, real yj,
857 real xk, real yk, real *a, real *b)
859 /* determine parameters by solving the equation system, since we know the
860 * virtual site coordinates here.
862 real dx_ij, dx_ik, dy_ij, dy_ik;
869 *a = ( (xd-xi)*dy_ik - dx_ik*(yd-yi) ) / (dx_ij*dy_ik - dx_ik*dy_ij);
870 *b = ( yd - yi - (*a)*dy_ij ) / dy_ik;
874 static int gen_vsites_trp(PreprocessingAtomTypes *atype,
875 std::vector<gmx::RVec> *newx,
876 t_atom *newatom[], char ***newatomname[],
877 int *o2n[], int *newvsite_type[], int *newcgnr[],
878 t_symtab *symtab, int *nadd,
879 gmx::ArrayRef<const gmx::RVec> x, int *cgnr[],
880 t_atoms *at, int *vsite_type[],
881 gmx::ArrayRef<InteractionTypeParameters> plist,
882 int nrfound, int *ats, int add_shift,
883 gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
886 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
888 atCB, atCG, atCD1, atHD1, atCD2, atNE1, atHE1, atCE2, atCE3, atHE3,
889 atCZ2, atHZ2, atCZ3, atHZ3, atCH2, atHH2, atNR
891 /* weights for determining the COM's of both rings (M1 and M2): */
892 real mw[NMASS][atNR] = {
893 { 0, 1, 1, 1, 0.5, 1, 1, 0.5, 0, 0,
895 { 0, 0, 0, 0, 0.5, 0, 0, 0.5, 1, 1,
899 real xi[atNR], yi[atNR];
900 real xcom[NMASS], ycom[NMASS], alpha;
901 real b_CD2_CE2, b_NE1_CE2, b_CG_CD2, b_CH2_HH2, b_CE2_CZ2;
902 real b_NE1_HE1, b_CD2_CE3, b_CE3_CZ3, b_CB_CG;
903 real b_CZ2_CH2, b_CZ2_HZ2, b_CD1_HD1, b_CE3_HE3;
904 real b_CG_CD1, b_CZ3_HZ3;
905 real a_NE1_CE2_CD2, a_CE2_CD2_CG, a_CB_CG_CD2, a_CE2_CD2_CE3;
906 real a_CD2_CG_CD1, a_CE2_CZ2_HZ2, a_CZ2_CH2_HH2;
907 real a_CD2_CE2_CZ2, a_CD2_CE3_CZ3, a_CE3_CZ3_HZ3, a_CG_CD1_HD1;
908 real a_CE2_CZ2_CH2, a_HE1_NE1_CE2, a_CD2_CE3_HE3;
909 int atM[NMASS], tpM, i, i0, j, nvsite;
910 real mM[NMASS], dCBM1, dCBM2, dM1M2;
912 rvec r_ij, r_ik, t1, t2;
917 gmx_incons("atom types in gen_vsites_trp");
919 /* Get geometry from database */
920 b_CD2_CE2 = get_ddb_bond(vsitetop, "TRP", "CD2", "CE2");
921 b_NE1_CE2 = get_ddb_bond(vsitetop, "TRP", "NE1", "CE2");
922 b_CG_CD1 = get_ddb_bond(vsitetop, "TRP", "CG", "CD1");
923 b_CG_CD2 = get_ddb_bond(vsitetop, "TRP", "CG", "CD2");
924 b_CB_CG = get_ddb_bond(vsitetop, "TRP", "CB", "CG");
925 b_CE2_CZ2 = get_ddb_bond(vsitetop, "TRP", "CE2", "CZ2");
926 b_CD2_CE3 = get_ddb_bond(vsitetop, "TRP", "CD2", "CE3");
927 b_CE3_CZ3 = get_ddb_bond(vsitetop, "TRP", "CE3", "CZ3");
928 b_CZ2_CH2 = get_ddb_bond(vsitetop, "TRP", "CZ2", "CH2");
930 b_CD1_HD1 = get_ddb_bond(vsitetop, "TRP", "CD1", "HD1");
931 b_CZ2_HZ2 = get_ddb_bond(vsitetop, "TRP", "CZ2", "HZ2");
932 b_NE1_HE1 = get_ddb_bond(vsitetop, "TRP", "NE1", "HE1");
933 b_CH2_HH2 = get_ddb_bond(vsitetop, "TRP", "CH2", "HH2");
934 b_CE3_HE3 = get_ddb_bond(vsitetop, "TRP", "CE3", "HE3");
935 b_CZ3_HZ3 = get_ddb_bond(vsitetop, "TRP", "CZ3", "HZ3");
937 a_NE1_CE2_CD2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "NE1", "CE2", "CD2");
938 a_CE2_CD2_CG = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CE2", "CD2", "CG");
939 a_CB_CG_CD2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CB", "CG", "CD2");
940 a_CD2_CG_CD1 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CD2", "CG", "CD1");
942 a_CE2_CD2_CE3 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CE2", "CD2", "CE3");
943 a_CD2_CE2_CZ2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CD2", "CE2", "CZ2");
944 a_CD2_CE3_CZ3 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CD2", "CE3", "CZ3");
945 a_CE3_CZ3_HZ3 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CE3", "CZ3", "HZ3");
946 a_CZ2_CH2_HH2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CZ2", "CH2", "HH2");
947 a_CE2_CZ2_HZ2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CE2", "CZ2", "HZ2");
948 a_CE2_CZ2_CH2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CE2", "CZ2", "CH2");
949 a_CG_CD1_HD1 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CG", "CD1", "HD1");
950 a_HE1_NE1_CE2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "HE1", "NE1", "CE2");
951 a_CD2_CE3_HE3 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CD2", "CE3", "HE3");
953 /* Calculate local coordinates.
954 * y-axis (x=0) is the bond CD2-CE2.
955 * x-axis (y=0) is perpendicular to the bond CD2-CE2 and
956 * intersects the middle of the bond.
959 yi[atCD2] = -0.5*b_CD2_CE2;
962 yi[atCE2] = 0.5*b_CD2_CE2;
964 xi[atNE1] = -b_NE1_CE2*std::sin(a_NE1_CE2_CD2);
965 yi[atNE1] = yi[atCE2]-b_NE1_CE2*std::cos(a_NE1_CE2_CD2);
967 xi[atCG] = -b_CG_CD2*std::sin(a_CE2_CD2_CG);
968 yi[atCG] = yi[atCD2]+b_CG_CD2*std::cos(a_CE2_CD2_CG);
970 alpha = a_CE2_CD2_CG + M_PI - a_CB_CG_CD2;
971 xi[atCB] = xi[atCG]-b_CB_CG*std::sin(alpha);
972 yi[atCB] = yi[atCG]+b_CB_CG*std::cos(alpha);
974 alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - M_PI;
975 xi[atCD1] = xi[atCG]-b_CG_CD1*std::sin(alpha);
976 yi[atCD1] = yi[atCG]+b_CG_CD1*std::cos(alpha);
978 xi[atCE3] = b_CD2_CE3*std::sin(a_CE2_CD2_CE3);
979 yi[atCE3] = yi[atCD2]+b_CD2_CE3*std::cos(a_CE2_CD2_CE3);
981 xi[atCZ2] = b_CE2_CZ2*std::sin(a_CD2_CE2_CZ2);
982 yi[atCZ2] = yi[atCE2]-b_CE2_CZ2*std::cos(a_CD2_CE2_CZ2);
984 alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - M_PI;
985 xi[atCZ3] = xi[atCE3]+b_CE3_CZ3*std::sin(alpha);
986 yi[atCZ3] = yi[atCE3]+b_CE3_CZ3*std::cos(alpha);
988 alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - M_PI;
989 xi[atCH2] = xi[atCZ2]+b_CZ2_CH2*std::sin(alpha);
990 yi[atCH2] = yi[atCZ2]-b_CZ2_CH2*std::cos(alpha);
993 alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - a_CG_CD1_HD1;
994 xi[atHD1] = xi[atCD1]-b_CD1_HD1*std::sin(alpha);
995 yi[atHD1] = yi[atCD1]+b_CD1_HD1*std::cos(alpha);
997 alpha = a_NE1_CE2_CD2 + M_PI - a_HE1_NE1_CE2;
998 xi[atHE1] = xi[atNE1]-b_NE1_HE1*std::sin(alpha);
999 yi[atHE1] = yi[atNE1]-b_NE1_HE1*std::cos(alpha);
1001 alpha = a_CE2_CD2_CE3 + M_PI - a_CD2_CE3_HE3;
1002 xi[atHE3] = xi[atCE3]+b_CE3_HE3*std::sin(alpha);
1003 yi[atHE3] = yi[atCE3]+b_CE3_HE3*std::cos(alpha);
1005 alpha = a_CD2_CE2_CZ2 + M_PI - a_CE2_CZ2_HZ2;
1006 xi[atHZ2] = xi[atCZ2]+b_CZ2_HZ2*std::sin(alpha);
1007 yi[atHZ2] = yi[atCZ2]-b_CZ2_HZ2*std::cos(alpha);
1009 alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - a_CZ2_CH2_HH2;
1010 xi[atHZ3] = xi[atCZ3]+b_CZ3_HZ3*std::sin(alpha);
1011 yi[atHZ3] = yi[atCZ3]+b_CZ3_HZ3*std::cos(alpha);
1013 alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - a_CE3_CZ3_HZ3;
1014 xi[atHH2] = xi[atCH2]+b_CH2_HH2*std::sin(alpha);
1015 yi[atHH2] = yi[atCH2]-b_CH2_HH2*std::cos(alpha);
1017 /* Calculate masses for each ring and put it on the dummy masses */
1018 for (j = 0; j < NMASS; j++)
1020 mM[j] = xcom[j] = ycom[j] = 0;
1022 for (i = 0; i < atNR; i++)
1026 for (j = 0; j < NMASS; j++)
1028 mM[j] += mw[j][i] * at->atom[ats[i]].m;
1029 xcom[j] += xi[i] * mw[j][i] * at->atom[ats[i]].m;
1030 ycom[j] += yi[i] * mw[j][i] * at->atom[ats[i]].m;
1034 for (j = 0; j < NMASS; j++)
1040 /* get dummy mass type */
1041 tpM = vsite_nm2type("MW", atype);
1042 /* make space for 2 masses: shift all atoms starting with CB */
1044 for (j = 0; j < NMASS; j++)
1046 atM[j] = i0+*nadd+j;
1050 fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, (*o2n)[i0]+1);
1053 for (j = i0; j < at->nr; j++)
1055 (*o2n)[j] = j+*nadd;
1057 newx->resize(at->nr+*nadd);
1058 srenew(*newatom, at->nr+*nadd);
1059 srenew(*newatomname, at->nr+*nadd);
1060 srenew(*newvsite_type, at->nr+*nadd);
1061 srenew(*newcgnr, at->nr+*nadd);
1062 for (j = 0; j < NMASS; j++)
1064 (*newatomname)[at->nr+*nadd-1-j] = nullptr;
1067 /* Dummy masses will be placed at the center-of-mass in each ring. */
1069 /* calc initial position for dummy masses in real (non-local) coordinates.
1070 * Cheat by using the routine to calculate virtual site parameters. It is
1071 * much easier when we have the coordinates expressed in terms of
1074 rvec_sub(x[ats[atCB]], x[ats[atCG]], r_ij);
1075 rvec_sub(x[ats[atCD2]], x[ats[atCG]], r_ik);
1076 calc_vsite3_param(xcom[0], ycom[0], xi[atCG], yi[atCG], xi[atCB], yi[atCB],
1077 xi[atCD2], yi[atCD2], &a, &b);
1080 rvec_add(t1, t2, t1);
1081 rvec_add(t1, x[ats[atCG]], (*newx)[atM[0]]);
1083 calc_vsite3_param(xcom[1], ycom[1], xi[atCG], yi[atCG], xi[atCB], yi[atCB],
1084 xi[atCD2], yi[atCD2], &a, &b);
1087 rvec_add(t1, t2, t1);
1088 rvec_add(t1, x[ats[atCG]], (*newx)[atM[1]]);
1090 /* set parameters for the masses */
1091 for (j = 0; j < NMASS; j++)
1093 sprintf(name, "MW%d", j+1);
1094 (*newatomname) [atM[j]] = put_symtab(symtab, name);
1095 (*newatom) [atM[j]].m = (*newatom)[atM[j]].mB = mM[j];
1096 (*newatom) [atM[j]].q = (*newatom)[atM[j]].qB = 0.0;
1097 (*newatom) [atM[j]].type = (*newatom)[atM[j]].typeB = tpM;
1098 (*newatom) [atM[j]].ptype = eptAtom;
1099 (*newatom) [atM[j]].resind = at->atom[i0].resind;
1100 (*newatom) [atM[j]].elem[0] = 'M';
1101 (*newatom) [atM[j]].elem[1] = '\0';
1102 (*newvsite_type)[atM[j]] = NOTSET;
1103 (*newcgnr) [atM[j]] = (*cgnr)[i0];
1105 /* renumber cgnr: */
1106 for (i = i0; i < at->nr; i++)
1111 /* constraints between CB, M1 and M2 */
1112 /* 'add_shift' says which atoms won't be renumbered afterwards */
1113 dCBM1 = std::hypot( xcom[0]-xi[atCB], ycom[0]-yi[atCB] );
1114 dM1M2 = std::hypot( xcom[0]-xcom[1], ycom[0]-ycom[1] );
1115 dCBM2 = std::hypot( xcom[1]-xi[atCB], ycom[1]-yi[atCB] );
1116 my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift+atM[0], dCBM1);
1117 my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift+atM[1], dCBM2);
1118 my_add_param(&(plist[F_CONSTRNC]), add_shift+atM[0], add_shift+atM[1], dM1M2);
1120 /* rest will be vsite3 */
1122 for (i = 0; i < atNR; i++)
1126 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
1127 (*vsite_type)[ats[i]] = F_VSITE3;
1132 /* now define all vsites from M1, M2, CB, ie:
1133 r_d = r_M1 + a r_M1_M2 + b r_M1_CB */
1134 for (i = 0; i < atNR; i++)
1136 if ( (*vsite_type)[ats[i]] == F_VSITE3)
1138 calc_vsite3_param(xi[i], yi[i], xcom[0], ycom[0], xcom[1], ycom[1], xi[atCB], yi[atCB], &a, &b);
1139 add_vsite3_param(&plist[F_VSITE3],
1140 ats[i], add_shift+atM[0], add_shift+atM[1], ats[atCB], a, b);
1148 static int gen_vsites_tyr(PreprocessingAtomTypes *atype,
1149 std::vector<gmx::RVec> *newx,
1150 t_atom *newatom[], char ***newatomname[],
1151 int *o2n[], int *newvsite_type[], int *newcgnr[],
1152 t_symtab *symtab, int *nadd,
1153 gmx::ArrayRef<const gmx::RVec> x, int *cgnr[],
1154 t_atoms *at, int *vsite_type[],
1155 gmx::ArrayRef<InteractionTypeParameters> plist,
1156 int nrfound, int *ats, int add_shift,
1157 gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
1159 int nvsite, i, i0, j, atM, tpM;
1160 real dCGCE, dCEOH, dCGM, tmp1, a, b;
1161 real bond_cc, bond_ch, bond_co, bond_oh, angle_coh;
1167 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1169 atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
1170 atCZ, atOH, atHH, atNR
1173 /* CG, CE1, CE2 (as in general 6-ring) and OH and HH stay,
1174 rest gets virtualized.
1175 Now we have two linked triangles with one improper keeping them flat */
1176 if (atNR != nrfound)
1178 gmx_incons("Number of atom types in gen_vsites_tyr");
1181 /* Aromatic rings have 6-fold symmetry, so we only need one bond length
1182 * for the ring part (angle is always 120 degrees).
1184 bond_cc = get_ddb_bond(vsitetop, "TYR", "CD1", "CE1");
1185 bond_ch = get_ddb_bond(vsitetop, "TYR", "CD1", "HD1");
1186 bond_co = get_ddb_bond(vsitetop, "TYR", "CZ", "OH");
1187 bond_oh = get_ddb_bond(vsitetop, "TYR", "OH", "HH");
1188 angle_coh = DEG2RAD*get_ddb_angle(vsitetop, "TYR", "CZ", "OH", "HH");
1190 xi[atCG] = -bond_cc+bond_cc*std::cos(ANGLE_6RING);
1191 xi[atCD1] = -bond_cc;
1192 xi[atHD1] = xi[atCD1]+bond_ch*std::cos(ANGLE_6RING);
1194 xi[atHE1] = xi[atCE1]-bond_ch*std::cos(ANGLE_6RING);
1195 xi[atCD2] = xi[atCD1];
1196 xi[atHD2] = xi[atHD1];
1197 xi[atCE2] = xi[atCE1];
1198 xi[atHE2] = xi[atHE1];
1199 xi[atCZ] = bond_cc*std::cos(0.5*ANGLE_6RING);
1200 xi[atOH] = xi[atCZ]+bond_co;
1203 for (i = 0; i < atOH; i++)
1205 xcom += xi[i]*at->atom[ats[i]].m;
1206 mtot += at->atom[ats[i]].m;
1210 /* first do 6 ring as default,
1211 except CZ (we'll do that different) and HZ (we don't have that): */
1212 nvsite = gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, FALSE);
1214 /* then construct CZ from the 2nd triangle */
1215 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
1216 a = b = 0.5 * bond_co / ( bond_co - bond_cc*std::cos(ANGLE_6RING) );
1217 add_vsite3_param(&plist[F_VSITE3],
1218 ats[atCZ], ats[atOH], ats[atCE1], ats[atCE2], a, b);
1219 at->atom[ats[atCZ]].m = at->atom[ats[atCZ]].mB = 0;
1221 /* constraints between CE1, CE2 and OH */
1222 dCGCE = std::sqrt( cosrule(bond_cc, bond_cc, ANGLE_6RING) );
1223 dCEOH = std::sqrt( cosrule(bond_cc, bond_co, ANGLE_6RING) );
1224 my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atOH], dCEOH);
1225 my_add_param(&(plist[F_CONSTRNC]), ats[atCE2], ats[atOH], dCEOH);
1227 /* We also want to constrain the angle C-O-H, but since CZ is constructed
1228 * we need to introduce a constraint to CG.
1229 * CG is much further away, so that will lead to instabilities in LINCS
1230 * when we constrain both CG-HH and OH-HH distances. Instead of requiring
1231 * the use of lincs_order=8 we introduce a dummy mass three times further
1232 * away from OH than HH. The mass is accordingly a third, with the remaining
1233 * 2/3 moved to OH. This shouldn't cause any problems since the forces will
1234 * apply to the HH constructed atom and not directly on the virtual mass.
1237 vdist = 2.0*bond_oh;
1238 mM = at->atom[ats[atHH]].m/2.0;
1239 at->atom[ats[atOH]].m += mM; /* add 1/2 of original H mass */
1240 at->atom[ats[atOH]].mB += mM; /* add 1/2 of original H mass */
1241 at->atom[ats[atHH]].m = at->atom[ats[atHH]].mB = 0;
1243 /* get dummy mass type */
1244 tpM = vsite_nm2type("MW", atype);
1245 /* make space for 1 mass: shift HH only */
1250 fprintf(stderr, "Inserting 1 dummy mass at %d\n", (*o2n)[i0]+1);
1253 for (j = i0; j < at->nr; j++)
1255 (*o2n)[j] = j+*nadd;
1257 newx->resize(at->nr+*nadd);
1258 srenew(*newatom, at->nr+*nadd);
1259 srenew(*newatomname, at->nr+*nadd);
1260 srenew(*newvsite_type, at->nr+*nadd);
1261 srenew(*newcgnr, at->nr+*nadd);
1262 (*newatomname)[at->nr+*nadd-1] = nullptr;
1264 /* Calc the dummy mass initial position */
1265 rvec_sub(x[ats[atHH]], x[ats[atOH]], r1);
1267 rvec_add(r1, x[ats[atHH]], (*newx)[atM]);
1269 strcpy(name, "MW1");
1270 (*newatomname) [atM] = put_symtab(symtab, name);
1271 (*newatom) [atM].m = (*newatom)[atM].mB = mM;
1272 (*newatom) [atM].q = (*newatom)[atM].qB = 0.0;
1273 (*newatom) [atM].type = (*newatom)[atM].typeB = tpM;
1274 (*newatom) [atM].ptype = eptAtom;
1275 (*newatom) [atM].resind = at->atom[i0].resind;
1276 (*newatom) [atM].elem[0] = 'M';
1277 (*newatom) [atM].elem[1] = '\0';
1278 (*newvsite_type)[atM] = NOTSET;
1279 (*newcgnr) [atM] = (*cgnr)[i0];
1280 /* renumber cgnr: */
1281 for (i = i0; i < at->nr; i++)
1286 (*vsite_type)[ats[atHH]] = F_VSITE2;
1288 /* assume we also want the COH angle constrained: */
1289 tmp1 = bond_cc*std::cos(0.5*ANGLE_6RING) + dCGCE*std::sin(ANGLE_6RING*0.5) + bond_co;
1290 dCGM = std::sqrt( cosrule(tmp1, vdist, angle_coh) );
1291 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], add_shift+atM, dCGM);
1292 my_add_param(&(plist[F_CONSTRNC]), ats[atOH], add_shift+atM, vdist);
1294 add_vsite2_param(&plist[F_VSITE2],
1295 ats[atHH], ats[atOH], add_shift+atM, 1.0/2.0);
1299 static int gen_vsites_his(t_atoms *at, int *vsite_type[],
1300 gmx::ArrayRef<InteractionTypeParameters> plist,
1301 int nrfound, int *ats, gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
1304 real a, b, alpha, dCGCE1, dCGNE2;
1305 real sinalpha, cosalpha;
1306 real xcom, ycom, mtot;
1307 real mG, mrest, mCE1, mNE2;
1308 real b_CG_ND1, b_ND1_CE1, b_CE1_NE2, b_CG_CD2, b_CD2_NE2;
1309 real b_ND1_HD1, b_NE2_HE2, b_CE1_HE1, b_CD2_HD2;
1310 real a_CG_ND1_CE1, a_CG_CD2_NE2, a_ND1_CE1_NE2, a_CE1_NE2_CD2;
1311 real a_NE2_CE1_HE1, a_NE2_CD2_HD2, a_CE1_ND1_HD1, a_CE1_NE2_HE2;
1314 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1316 atCG, atND1, atHD1, atCD2, atHD2, atCE1, atHE1, atNE2, atHE2, atNR
1318 real x[atNR], y[atNR];
1320 /* CG, CE1 and NE2 stay, each gets part of the total mass,
1321 rest gets virtualized */
1322 /* check number of atoms, 3 hydrogens may be missing: */
1323 /* assert( nrfound >= atNR-3 || nrfound <= atNR );
1324 * Don't understand the above logic. Shouldn't it be && rather than || ???
1326 if ((nrfound < atNR-3) || (nrfound > atNR))
1328 gmx_incons("Generating vsites for HIS");
1331 /* avoid warnings about uninitialized variables */
1332 b_ND1_HD1 = b_NE2_HE2 = b_CE1_HE1 = b_CD2_HD2 = a_NE2_CE1_HE1 =
1333 a_NE2_CD2_HD2 = a_CE1_ND1_HD1 = a_CE1_NE2_HE2 = 0;
1335 if (ats[atHD1] != NOTSET)
1337 if (ats[atHE2] != NOTSET)
1339 sprintf(resname, "HISH");
1343 sprintf(resname, "HISA");
1348 sprintf(resname, "HISB");
1351 /* Get geometry from database */
1352 b_CG_ND1 = get_ddb_bond(vsitetop, resname, "CG", "ND1");
1353 b_ND1_CE1 = get_ddb_bond(vsitetop, resname, "ND1", "CE1");
1354 b_CE1_NE2 = get_ddb_bond(vsitetop, resname, "CE1", "NE2");
1355 b_CG_CD2 = get_ddb_bond(vsitetop, resname, "CG", "CD2");
1356 b_CD2_NE2 = get_ddb_bond(vsitetop, resname, "CD2", "NE2");
1357 a_CG_ND1_CE1 = DEG2RAD*get_ddb_angle(vsitetop, resname, "CG", "ND1", "CE1");
1358 a_CG_CD2_NE2 = DEG2RAD*get_ddb_angle(vsitetop, resname, "CG", "CD2", "NE2");
1359 a_ND1_CE1_NE2 = DEG2RAD*get_ddb_angle(vsitetop, resname, "ND1", "CE1", "NE2");
1360 a_CE1_NE2_CD2 = DEG2RAD*get_ddb_angle(vsitetop, resname, "CE1", "NE2", "CD2");
1362 if (ats[atHD1] != NOTSET)
1364 b_ND1_HD1 = get_ddb_bond(vsitetop, resname, "ND1", "HD1");
1365 a_CE1_ND1_HD1 = DEG2RAD*get_ddb_angle(vsitetop, resname, "CE1", "ND1", "HD1");
1367 if (ats[atHE2] != NOTSET)
1369 b_NE2_HE2 = get_ddb_bond(vsitetop, resname, "NE2", "HE2");
1370 a_CE1_NE2_HE2 = DEG2RAD*get_ddb_angle(vsitetop, resname, "CE1", "NE2", "HE2");
1372 if (ats[atHD2] != NOTSET)
1374 b_CD2_HD2 = get_ddb_bond(vsitetop, resname, "CD2", "HD2");
1375 a_NE2_CD2_HD2 = DEG2RAD*get_ddb_angle(vsitetop, resname, "NE2", "CD2", "HD2");
1377 if (ats[atHE1] != NOTSET)
1379 b_CE1_HE1 = get_ddb_bond(vsitetop, resname, "CE1", "HE1");
1380 a_NE2_CE1_HE1 = DEG2RAD*get_ddb_angle(vsitetop, resname, "NE2", "CE1", "HE1");
1383 /* constraints between CG, CE1 and NE1 */
1384 dCGCE1 = std::sqrt( cosrule(b_CG_ND1, b_ND1_CE1, a_CG_ND1_CE1) );
1385 dCGNE2 = std::sqrt( cosrule(b_CG_CD2, b_CD2_NE2, a_CG_CD2_NE2) );
1387 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE1);
1388 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atNE2], dCGNE2);
1389 /* we already have a constraint CE1-NE2, so we don't add it again */
1391 /* calculate the positions in a local frame of reference.
1392 * The x-axis is the line from CG that makes a right angle
1393 * with the bond CE1-NE2, and the y-axis the bond CE1-NE2.
1395 /* First calculate the x-axis intersection with y-axis (=yCE1).
1396 * Get cos(angle CG-CE1-NE2) :
1398 cosalpha = acosrule(dCGNE2, dCGCE1, b_CE1_NE2);
1400 y[atCE1] = cosalpha*dCGCE1;
1402 y[atNE2] = y[atCE1]-b_CE1_NE2;
1403 sinalpha = std::sqrt(1-cosalpha*cosalpha);
1404 x[atCG] = -sinalpha*dCGCE1;
1406 x[atHE1] = x[atHE2] = x[atHD1] = x[atHD2] = 0;
1407 y[atHE1] = y[atHE2] = y[atHD1] = y[atHD2] = 0;
1409 /* calculate ND1 and CD2 positions from CE1 and NE2 */
1411 x[atND1] = -b_ND1_CE1*std::sin(a_ND1_CE1_NE2);
1412 y[atND1] = y[atCE1]-b_ND1_CE1*std::cos(a_ND1_CE1_NE2);
1414 x[atCD2] = -b_CD2_NE2*std::sin(a_CE1_NE2_CD2);
1415 y[atCD2] = y[atNE2]+b_CD2_NE2*std::cos(a_CE1_NE2_CD2);
1417 /* And finally the hydrogen positions */
1418 if (ats[atHE1] != NOTSET)
1420 x[atHE1] = x[atCE1] + b_CE1_HE1*std::sin(a_NE2_CE1_HE1);
1421 y[atHE1] = y[atCE1] - b_CE1_HE1*std::cos(a_NE2_CE1_HE1);
1423 /* HD2 - first get (ccw) angle from (positive) y-axis */
1424 if (ats[atHD2] != NOTSET)
1426 alpha = a_CE1_NE2_CD2 + M_PI - a_NE2_CD2_HD2;
1427 x[atHD2] = x[atCD2] - b_CD2_HD2*std::sin(alpha);
1428 y[atHD2] = y[atCD2] + b_CD2_HD2*std::cos(alpha);
1430 if (ats[atHD1] != NOTSET)
1432 /* HD1 - first get (cw) angle from (positive) y-axis */
1433 alpha = a_ND1_CE1_NE2 + M_PI - a_CE1_ND1_HD1;
1434 x[atHD1] = x[atND1] - b_ND1_HD1*std::sin(alpha);
1435 y[atHD1] = y[atND1] - b_ND1_HD1*std::cos(alpha);
1437 if (ats[atHE2] != NOTSET)
1439 x[atHE2] = x[atNE2] + b_NE2_HE2*std::sin(a_CE1_NE2_HE2);
1440 y[atHE2] = y[atNE2] + b_NE2_HE2*std::cos(a_CE1_NE2_HE2);
1442 /* Have all coordinates now */
1444 /* calc center-of-mass; keep atoms CG, CE1, NE2 and
1445 * set the rest to vsite3
1447 mtot = xcom = ycom = 0;
1449 for (i = 0; i < atNR; i++)
1451 if (ats[i] != NOTSET)
1453 mtot += at->atom[ats[i]].m;
1454 xcom += x[i]*at->atom[ats[i]].m;
1455 ycom += y[i]*at->atom[ats[i]].m;
1456 if (i != atCG && i != atCE1 && i != atNE2)
1458 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
1459 (*vsite_type)[ats[i]] = F_VSITE3;
1464 if (nvsite+3 != nrfound)
1466 gmx_incons("Generating vsites for HIS");
1472 /* distribute mass so that com stays the same */
1473 mG = xcom*mtot/x[atCG];
1475 mCE1 = (ycom-y[atNE2])*mrest/(y[atCE1]-y[atNE2]);
1478 at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = mG;
1479 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB = mCE1;
1480 at->atom[ats[atNE2]].m = at->atom[ats[atNE2]].mB = mNE2;
1483 if (ats[atHE1] != NOTSET)
1485 calc_vsite3_param(x[atHE1], y[atHE1], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1486 x[atCG], y[atCG], &a, &b);
1487 add_vsite3_param(&plist[F_VSITE3],
1488 ats[atHE1], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1491 if (ats[atHE2] != NOTSET)
1493 calc_vsite3_param(x[atHE2], y[atHE2], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1494 x[atCG], y[atCG], &a, &b);
1495 add_vsite3_param(&plist[F_VSITE3],
1496 ats[atHE2], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1500 calc_vsite3_param(x[atND1], y[atND1], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1501 x[atCG], y[atCG], &a, &b);
1502 add_vsite3_param(&plist[F_VSITE3],
1503 ats[atND1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1506 calc_vsite3_param(x[atCD2], y[atCD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1507 x[atCG], y[atCG], &a, &b);
1508 add_vsite3_param(&plist[F_VSITE3],
1509 ats[atCD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1512 if (ats[atHD1] != NOTSET)
1514 calc_vsite3_param(x[atHD1], y[atHD1], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1515 x[atCG], y[atCG], &a, &b);
1516 add_vsite3_param(&plist[F_VSITE3],
1517 ats[atHD1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1520 if (ats[atHD2] != NOTSET)
1522 calc_vsite3_param(x[atHD2], y[atHD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1523 x[atCG], y[atCG], &a, &b);
1524 add_vsite3_param(&plist[F_VSITE3],
1525 ats[atHD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1530 static bool is_vsite(int vsite_type)
1532 if (vsite_type == NOTSET)
1536 switch (abs(vsite_type) )
1550 static char atomnamesuffix[] = "1234";
1552 void do_vsites(gmx::ArrayRef<const PreprocessResidue> rtpFFDB, PreprocessingAtomTypes *atype,
1553 t_atoms *at, t_symtab *symtab,
1554 std::vector<gmx::RVec> *x,
1555 gmx::ArrayRef<InteractionTypeParameters> plist, int *vsite_type[], int *cgnr[],
1556 real mHmult, bool bVsiteAromatics,
1559 #define MAXATOMSPERRESIDUE 16
1560 int k, m, i0, ni0, whatres, resind, add_shift, nvsite, nadd;
1562 int nrfound = 0, needed, nrbonds, nrHatoms, Heavy, nrheavies, tpM, tpHeavy;
1563 int Hatoms[4], heavies[4];
1564 bool bWARNING, bAddVsiteParam, bFirstWater;
1566 real mHtot, mtot, fact, fact2;
1567 rvec rpar, rperp, temp;
1568 char tpname[32], nexttpname[32];
1569 int *o2n, *newvsite_type, *newcgnr, ats[MAXATOMSPERRESIDUE];
1571 InteractionTypeParameters *params;
1572 char ***newatomname;
1573 char *resnm = nullptr;
1575 bool isN, planarN, bFound;
1577 /* if bVsiteAromatics=TRUE do_vsites will specifically convert atoms in
1578 PHE, TRP, TYR and HIS to a construction of virtual sites */
1580 resPHE, resTRP, resTYR, resHIS, resNR
1582 const char *resnms[resNR] = { "PHE", "TRP", "TYR", "HIS" };
1583 /* Amber03 alternative names for termini */
1584 const char *resnmsN[resNR] = { "NPHE", "NTRP", "NTYR", "NHIS" };
1585 const char *resnmsC[resNR] = { "CPHE", "CTRP", "CTYR", "CHIS" };
1586 /* HIS can be known as HISH, HIS1, HISA, HID, HIE, HIP, etc. too */
1587 bool bPartial[resNR] = { FALSE, FALSE, FALSE, TRUE };
1588 /* the atnms for every residue MUST correspond to the enums in the
1589 gen_vsites_* (one for each residue) routines! */
1590 /* also the atom names in atnms MUST be in the same order as in the .rtp! */
1591 const char *atnms[resNR][MAXATOMSPERRESIDUE+1] = {
1593 "CD1", "HD1", "CD2", "HD2",
1594 "CE1", "HE1", "CE2", "HE2",
1595 "CZ", "HZ", nullptr },
1598 "CD1", "HD1", "CD2",
1599 "NE1", "HE1", "CE2", "CE3", "HE3",
1600 "CZ2", "HZ2", "CZ3", "HZ3",
1601 "CH2", "HH2", nullptr },
1603 "CD1", "HD1", "CD2", "HD2",
1604 "CE1", "HE1", "CE2", "HE2",
1605 "CZ", "OH", "HH", nullptr },
1607 "ND1", "HD1", "CD2", "HD2",
1608 "CE1", "HE1", "NE2", "HE2", nullptr }
1613 printf("Searching for atoms to make virtual sites ...\n");
1614 fprintf(debug, "# # # VSITES # # #\n");
1617 std::vector<std::string> db = fflib_search_file_end(ffdir, ".vsd", FALSE);
1619 /* Container of CH3/NH3/NH2 configuration entries.
1620 * See comments in read_vsite_database. It isnt beautiful,
1621 * but it had to be fixed, and I dont even want to try to
1622 * maintain this part of the code...
1624 std::vector<VirtualSiteConfiguration> vsiteconflist;
1626 // TODO those have been deprecated and should be removed completely.
1627 /* Container of geometry (bond/angle) entries for
1628 * residues like PHE, TRP, TYR, HIS, etc., where we need
1629 * to know the geometry to construct vsite aromatics.
1630 * Note that equilibrium geometry isnt necessarily the same
1631 * as the individual bond and angle values given in the
1632 * force field (rings can be strained).
1634 std::vector<VirtualSiteTopology> vsitetop;
1635 for (const auto &filename : db)
1637 read_vsite_database(filename.c_str(), &vsiteconflist, &vsitetop);
1643 /* we need a marker for which atoms should *not* be renumbered afterwards */
1644 add_shift = 10*at->nr;
1645 /* make arrays where masses can be inserted into */
1646 std::vector<gmx::RVec> newx(at->nr);
1647 snew(newatom, at->nr);
1648 snew(newatomname, at->nr);
1649 snew(newvsite_type, at->nr);
1650 snew(newcgnr, at->nr);
1651 /* make index array to tell where the atoms go to when masses are inserted */
1653 for (int i = 0; i < at->nr; i++)
1657 /* make index to tell which residues were already processed */
1658 std::vector<bool> bResProcessed(at->nres);
1662 /* generate vsite constructions */
1663 /* loop over all atoms */
1665 for (int i = 0; (i < at->nr); i++)
1667 if (at->atom[i].resind != resind)
1669 resind = at->atom[i].resind;
1670 resnm = *(at->resinfo[resind].name);
1672 /* first check for aromatics to virtualize */
1673 /* don't waste our effort on DNA, water etc. */
1674 /* Only do the vsite aromatic stuff when we reach the
1675 * CA atom, since there might be an X2/X3 group on the
1676 * N-terminus that must be treated first.
1678 if (bVsiteAromatics &&
1679 (strcmp(*(at->atomname[i]), "CA") == 0) &&
1680 !bResProcessed[resind] &&
1681 rt.namedResidueHasType(*(at->resinfo[resind].name), "Protein") )
1683 /* mark this residue */
1684 bResProcessed[resind] = TRUE;
1685 /* find out if this residue needs converting */
1687 for (int j = 0; j < resNR && whatres == NOTSET; j++)
1690 cmplength = bPartial[j] ? strlen(resnm)-1 : strlen(resnm);
1692 bFound = ((gmx_strncasecmp(resnm, resnms[j], cmplength) == 0) ||
1693 (gmx_strncasecmp(resnm, resnmsN[j], cmplength) == 0) ||
1694 (gmx_strncasecmp(resnm, resnmsC[j], cmplength) == 0));
1699 /* get atoms we will be needing for the conversion */
1701 for (k = 0; atnms[j][k]; k++)
1704 for (m = i; m < at->nr && at->atom[m].resind == resind && ats[k] == NOTSET; m++)
1706 if (gmx_strcasecmp(*(at->atomname[m]), atnms[j][k]) == 0)
1714 /* now k is number of atom names in atnms[j] */
1723 if (nrfound < needed)
1725 gmx_fatal(FARGS, "not enough atoms found (%d, need %d) in "
1726 "residue %s %d while\n "
1727 "generating aromatics virtual site construction",
1728 nrfound, needed, resnm, at->resinfo[resind].nr);
1730 /* Advance overall atom counter */
1734 /* the enums for every residue MUST correspond to atnms[residue] */
1740 fprintf(stderr, "PHE at %d\n", o2n[ats[0]]+1);
1742 nvsite += gen_vsites_phe(at, vsite_type, plist, nrfound, ats, vsitetop);
1747 fprintf(stderr, "TRP at %d\n", o2n[ats[0]]+1);
1749 nvsite += gen_vsites_trp(atype, &newx, &newatom, &newatomname, &o2n,
1750 &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
1751 at, vsite_type, plist, nrfound, ats, add_shift, vsitetop);
1756 fprintf(stderr, "TYR at %d\n", o2n[ats[0]]+1);
1758 nvsite += gen_vsites_tyr(atype, &newx, &newatom, &newatomname, &o2n,
1759 &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
1760 at, vsite_type, plist, nrfound, ats, add_shift, vsitetop);
1765 fprintf(stderr, "HIS at %d\n", o2n[ats[0]]+1);
1767 nvsite += gen_vsites_his(at, vsite_type, plist, nrfound, ats, vsitetop);
1770 /* this means this residue won't be processed */
1773 gmx_fatal(FARGS, "DEATH HORROR in do_vsites (%s:%d)",
1774 __FILE__, __LINE__);
1775 } /* switch whatres */
1776 /* skip back to beginning of residue */
1777 while (i > 0 && at->atom[i-1].resind == resind)
1781 } /* if bVsiteAromatics & is protein */
1783 /* now process the rest of the hydrogens */
1784 /* only process hydrogen atoms which are not already set */
1785 if ( ((*vsite_type)[i] == NOTSET) && is_hydrogen(*(at->atomname[i])))
1787 /* find heavy atom, count #bonds from it and #H atoms bound to it
1788 and return H atom numbers (Hatoms) and heavy atom numbers (heavies) */
1789 count_bonds(i, &plist[F_BONDS], at->atomname,
1790 &nrbonds, &nrHatoms, Hatoms, &Heavy, &nrheavies, heavies);
1791 /* get Heavy atom type */
1792 tpHeavy = get_atype(Heavy, at, rtpFFDB, &rt);
1793 strcpy(tpname, atype->atomNameFromAtomType(tpHeavy));
1796 bAddVsiteParam = TRUE;
1797 /* nested if's which check nrHatoms, nrbonds and atomname */
1803 (*vsite_type)[i] = F_BONDS;
1805 case 3: /* =CH-, -NH- or =NH+- */
1806 (*vsite_type)[i] = F_VSITE3FD;
1808 case 4: /* --CH- (tert) */
1809 /* The old type 4FD had stability issues, so
1810 * all new constructs should use 4FDN
1812 (*vsite_type)[i] = F_VSITE4FDN;
1814 /* Check parity of heavy atoms from coordinates */
1819 rvec_sub((*x)[aj], (*x)[ai], tmpmat[0]);
1820 rvec_sub((*x)[ak], (*x)[ai], tmpmat[1]);
1821 rvec_sub((*x)[al], (*x)[ai], tmpmat[2]);
1823 if (det(tmpmat) > 0)
1831 default: /* nrbonds != 2, 3 or 4 */
1836 else if ( (nrHatoms == 2) && (nrbonds == 2) &&
1837 (at->atom[Heavy].atomnumber == 8) )
1839 bAddVsiteParam = FALSE; /* this is water: skip these hydrogens */
1842 bFirstWater = FALSE;
1846 "Not converting hydrogens in water to virtual sites\n");
1850 else if ( (nrHatoms == 2) && (nrbonds == 4) )
1852 /* -CH2- , -NH2+- */
1853 (*vsite_type)[Hatoms[0]] = F_VSITE3OUT;
1854 (*vsite_type)[Hatoms[1]] = -F_VSITE3OUT;
1858 /* 2 or 3 hydrogen atom, with 3 or 4 bonds in total to the heavy atom.
1859 * If it is a nitrogen, first check if it is planar.
1861 isN = planarN = FALSE;
1862 if ((nrHatoms == 2) && ((*at->atomname[Heavy])[0] == 'N'))
1865 int j = nitrogen_is_planar(vsiteconflist, tpname);
1868 gmx_fatal(FARGS, "No vsite database NH2 entry for type %s\n", tpname);
1872 if ( (nrHatoms == 2) && (nrbonds == 3) && ( !isN || planarN ) )
1874 /* =CH2 or, if it is a nitrogen NH2, it is a planar one */
1875 (*vsite_type)[Hatoms[0]] = F_VSITE3FAD;
1876 (*vsite_type)[Hatoms[1]] = -F_VSITE3FAD;
1878 else if ( ( (nrHatoms == 2) && (nrbonds == 3) &&
1879 ( isN && !planarN ) ) ||
1880 ( (nrHatoms == 3) && (nrbonds == 4) ) )
1882 /* CH3, NH3 or non-planar NH2 group */
1883 int Hat_vsite_type[3] = { F_VSITE3, F_VSITE3OUT, F_VSITE3OUT };
1884 bool Hat_SwapParity[3] = { FALSE, TRUE, FALSE };
1888 fprintf(stderr, "-XH3 or nonplanar NH2 group at %d\n", i+1);
1890 bAddVsiteParam = FALSE; /* we'll do this ourselves! */
1891 /* -NH2 (umbrella), -NH3+ or -CH3 */
1892 (*vsite_type)[Heavy] = F_VSITE3;
1893 for (int j = 0; j < nrHatoms; j++)
1895 (*vsite_type)[Hatoms[j]] = Hat_vsite_type[j];
1897 /* get dummy mass type from first char of heavy atom type (N or C) */
1899 strcpy(nexttpname, atype->atomNameFromAtomType(get_atype(heavies[0], at, rtpFFDB, &rt)));
1900 std::string ch = get_dummymass_name(vsiteconflist, tpname, nexttpname);
1906 gmx_fatal(FARGS, "Can't find dummy mass for type %s bonded to type %s in the virtual site database (.vsd files). Add it to the database!\n", tpname, nexttpname);
1910 gmx_fatal(FARGS, "A dummy mass for type %s bonded to type %s is required, but no virtual site database (.vsd) files where found.\n", tpname, nexttpname);
1918 tpM = vsite_nm2type(name.c_str(), atype);
1919 /* make space for 2 masses: shift all atoms starting with 'Heavy' */
1925 fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, o2n[i0]+1);
1928 for (int j = i0; j < at->nr; j++)
1933 newx.resize(at->nr+nadd);
1934 srenew(newatom, at->nr+nadd);
1935 srenew(newatomname, at->nr+nadd);
1936 srenew(newvsite_type, at->nr+nadd);
1937 srenew(newcgnr, at->nr+nadd);
1939 for (int j = 0; j < NMASS; j++)
1941 newatomname[at->nr+nadd-1-j] = nullptr;
1944 /* calculate starting position for the masses */
1946 /* get atom masses, and set Heavy and Hatoms mass to zero */
1947 for (int j = 0; j < nrHatoms; j++)
1949 mHtot += get_amass(Hatoms[j], at, rtpFFDB, &rt);
1950 at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
1952 mtot = mHtot + get_amass(Heavy, at, rtpFFDB, &rt);
1953 at->atom[Heavy].m = at->atom[Heavy].mB = 0;
1959 fact = std::sqrt(fact2);
1960 /* generate vectors parallel and perpendicular to rotational axis:
1961 * rpar = Heavy -> Hcom
1962 * rperp = Hcom -> H1 */
1964 for (int j = 0; j < nrHatoms; j++)
1966 rvec_inc(rpar, (*x)[Hatoms[j]]);
1968 svmul(1.0/nrHatoms, rpar, rpar); /* rpar = ( H1+H2+H3 ) / 3 */
1969 rvec_dec(rpar, (*x)[Heavy]); /* - Heavy */
1970 rvec_sub((*x)[Hatoms[0]], (*x)[Heavy], rperp);
1971 rvec_dec(rperp, rpar); /* rperp = H1 - Heavy - rpar */
1972 /* calc mass positions */
1973 svmul(fact2, rpar, temp);
1974 for (int j = 0; (j < NMASS); j++) /* xM = xN + fact2 * rpar +/- fact * rperp */
1976 rvec_add((*x)[Heavy], temp, newx[ni0+j]);
1978 svmul(fact, rperp, temp);
1979 rvec_inc(newx[ni0 ], temp);
1980 rvec_dec(newx[ni0+1], temp);
1981 /* set atom parameters for the masses */
1982 for (int j = 0; (j < NMASS); j++)
1984 /* make name: "M??#" or "M?#" (? is atomname, # is number) */
1987 for (k = 0; (*at->atomname[Heavy])[k] && ( k < NMASS ); k++)
1989 name[k+1] = (*at->atomname[Heavy])[k];
1991 name[k+1] = atomnamesuffix[j];
1993 newatomname[ni0+j] = put_symtab(symtab, name.c_str());
1994 newatom[ni0+j].m = newatom[ni0+j].mB = mtot/NMASS;
1995 newatom[ni0+j].q = newatom[ni0+j].qB = 0.0;
1996 newatom[ni0+j].type = newatom[ni0+j].typeB = tpM;
1997 newatom[ni0+j].ptype = eptAtom;
1998 newatom[ni0+j].resind = at->atom[i0].resind;
1999 newatom[ni0+j].elem[0] = 'M';
2000 newatom[ni0+j].elem[1] = '\0';
2001 newvsite_type[ni0+j] = NOTSET;
2002 newcgnr[ni0+j] = (*cgnr)[i0];
2004 /* add constraints between dummy masses and to heavies[0] */
2005 /* 'add_shift' says which atoms won't be renumbered afterwards */
2006 my_add_param(&(plist[F_CONSTRNC]), heavies[0], add_shift+ni0, NOTSET);
2007 my_add_param(&(plist[F_CONSTRNC]), heavies[0], add_shift+ni0+1, NOTSET);
2008 my_add_param(&(plist[F_CONSTRNC]), add_shift+ni0, add_shift+ni0+1, NOTSET);
2010 /* generate Heavy, H1, H2 and H3 from M1, M2 and heavies[0] */
2011 /* note that vsite_type cannot be NOTSET, because we just set it */
2012 add_vsite3_atoms (&plist[(*vsite_type)[Heavy]],
2013 Heavy, heavies[0], add_shift+ni0, add_shift+ni0+1,
2015 for (int j = 0; j < nrHatoms; j++)
2017 add_vsite3_atoms(&plist[(*vsite_type)[Hatoms[j]]],
2018 Hatoms[j], heavies[0], add_shift+ni0, add_shift+ni0+1,
2032 "Warning: cannot convert atom %d %s (bound to a heavy atom "
2034 " %d bonds and %d bound hydrogens atoms) to virtual site\n",
2035 i+1, *(at->atomname[i]), tpname, nrbonds, nrHatoms);
2039 /* add vsite parameters to topology,
2040 also get rid of negative vsite_types */
2041 add_vsites(plist, (*vsite_type), Heavy, nrHatoms, Hatoms,
2042 nrheavies, heavies);
2043 /* transfer mass of virtual site to Heavy atom */
2044 for (int j = 0; j < nrHatoms; j++)
2046 if (is_vsite((*vsite_type)[Hatoms[j]]))
2048 at->atom[Heavy].m += at->atom[Hatoms[j]].m;
2049 at->atom[Heavy].mB = at->atom[Heavy].m;
2050 at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
2057 fprintf(debug, "atom %d: ", o2n[i]+1);
2058 print_bonds(debug, o2n, nrHatoms, Hatoms, Heavy, nrheavies, heavies);
2060 } /* if vsite NOTSET & is hydrogen */
2062 } /* for i < at->nr */
2066 fprintf(debug, "Before inserting new atoms:\n");
2067 for (int i = 0; i < at->nr; i++)
2069 fprintf(debug, "%4d %4d %4s %4d %4s %6d %-10s\n", i+1, o2n[i]+1,
2070 at->atomname[i] ? *(at->atomname[i]) : "(NULL)",
2071 at->resinfo[at->atom[i].resind].nr,
2072 at->resinfo[at->atom[i].resind].name ?
2073 *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
2075 ((*vsite_type)[i] == NOTSET) ?
2076 "NOTSET" : interaction_function[(*vsite_type)[i]].name);
2078 fprintf(debug, "new atoms to be inserted:\n");
2079 for (int i = 0; i < at->nr+nadd; i++)
2083 fprintf(debug, "%4d %4s %4d %6d %-10s\n", i+1,
2084 newatomname[i] ? *(newatomname[i]) : "(NULL)",
2085 newatom[i].resind, newcgnr[i],
2086 (newvsite_type[i] == NOTSET) ?
2087 "NOTSET" : interaction_function[newvsite_type[i]].name);
2092 /* add all original atoms to the new arrays, using o2n index array */
2093 for (int i = 0; i < at->nr; i++)
2095 newatomname [o2n[i]] = at->atomname [i];
2096 newatom [o2n[i]] = at->atom [i];
2097 newvsite_type[o2n[i]] = (*vsite_type)[i];
2098 newcgnr [o2n[i]] = (*cgnr) [i];
2099 copy_rvec((*x)[i], newx[o2n[i]]);
2101 /* throw away old atoms */
2103 sfree(at->atomname);
2106 /* put in the new ones */
2109 at->atomname = newatomname;
2110 *vsite_type = newvsite_type;
2113 if (at->nr > add_shift)
2115 gmx_fatal(FARGS, "Added impossible amount of dummy masses "
2116 "(%d on a total of %d atoms)\n", nadd, at->nr-nadd);
2121 fprintf(debug, "After inserting new atoms:\n");
2122 for (int i = 0; i < at->nr; i++)
2124 fprintf(debug, "%4d %4s %4d %4s %6d %-10s\n", i+1,
2125 at->atomname[i] ? *(at->atomname[i]) : "(NULL)",
2126 at->resinfo[at->atom[i].resind].nr,
2127 at->resinfo[at->atom[i].resind].name ?
2128 *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
2130 ((*vsite_type)[i] == NOTSET) ?
2131 "NOTSET" : interaction_function[(*vsite_type)[i]].name);
2135 /* now renumber all the interactions because of the added atoms */
2136 for (int ftype = 0; ftype < F_NRE; ftype++)
2138 InteractionTypeParameters *params = &(plist[ftype]);
2141 fprintf(debug, "Renumbering %d %s\n", params->nr,
2142 interaction_function[ftype].longname);
2144 for (int i = 0; i < params->nr; i++)
2146 for (int j = 0; j < NRAL(ftype); j++)
2148 if (params->param[i].a[j] >= add_shift)
2152 fprintf(debug, " [%d -> %d]", params->param[i].a[j],
2153 params->param[i].a[j]-add_shift);
2155 params->param[i].a[j] = params->param[i].a[j]-add_shift;
2161 fprintf(debug, " [%d -> %d]", params->param[i].a[j],
2162 o2n[params->param[i].a[j]]);
2164 params->param[i].a[j] = o2n[params->param[i].a[j]];
2169 fprintf(debug, "\n");
2173 /* now check if atoms in the added constraints are in increasing order */
2174 params = &(plist[F_CONSTRNC]);
2175 for (int i = 0; i < params->nr; i++)
2177 if (params->param[i].ai() > params->param[i].aj())
2179 int j = params->param[i].aj();
2180 params->param[i].aj() = params->param[i].ai();
2181 params->param[i].ai() = j;
2188 /* tell the user what we did */
2189 fprintf(stderr, "Marked %d virtual sites\n", nvsite);
2190 fprintf(stderr, "Added %d dummy masses\n", nadd);
2191 fprintf(stderr, "Added %d new constraints\n", plist[F_CONSTRNC].nr);
2194 void do_h_mass(InteractionTypeParameters *psb, int vsite_type[], t_atoms *at, real mHmult,
2197 /* loop over all atoms */
2198 for (int i = 0; i < at->nr; i++)
2200 /* adjust masses if i is hydrogen and not a virtual site */
2201 if (!is_vsite(vsite_type[i]) && is_hydrogen(*(at->atomname[i])) )
2203 /* find bonded heavy atom */
2205 for (int j = 0; (j < psb->nr) && (a == NOTSET); j++)
2207 /* if other atom is not a virtual site, it is the one we want */
2208 if ( (psb->param[j].ai() == i) &&
2209 !is_vsite(vsite_type[psb->param[j].aj()]) )
2211 a = psb->param[j].aj();
2213 else if ( (psb->param[j].aj() == i) &&
2214 !is_vsite(vsite_type[psb->param[j].ai()]) )
2216 a = psb->param[j].ai();
2221 gmx_fatal(FARGS, "Unbound hydrogen atom (%d) found while adjusting mass",
2225 /* adjust mass of i (hydrogen) with mHmult
2226 and correct mass of a (bonded atom) with same amount */
2229 at->atom[a].m -= (mHmult-1.0)*at->atom[i].m;
2230 at->atom[a].mB -= (mHmult-1.0)*at->atom[i].m;
2232 at->atom[i].m *= mHmult;
2233 at->atom[i].mB *= mHmult;