#include "resall.h"
#define MAXNAME 32
-#define OPENDIR '[' /* starting sign for directive */
-#define CLOSEDIR ']' /* ending sign for directive */
+#define OPENDIR '[' /* starting sign for directive */
+#define CLOSEDIR ']' /* ending sign for directive */
/*! \libinternal \brief
* The configuration describing a virtual site.
* \param[in] nextheavy Type of bonded heavy atom.
* \param[in] dummy What kind of dummy is used in the vsite.
*/
- explicit VirtualSiteConfiguration(const std::string &type, bool planar,
- int nhyd, const std::string &nextheavy, const std::string &dummy)
- : atomtype(type), isplanar(planar), nHydrogens(nhyd), nextHeavyType(nextheavy),
- dummyMass(dummy)
- {}
+ explicit VirtualSiteConfiguration(const std::string& type,
+ bool planar,
+ int nhyd,
+ const std::string& nextheavy,
+ const std::string& dummy) :
+ atomtype(type),
+ isplanar(planar),
+ nHydrogens(nhyd),
+ nextHeavyType(nextheavy),
+ dummyMass(dummy)
+ {
+ }
//! Type for the XH3/XH2 atom.
std::string atomtype;
/*! \brief Is the configuration planar?
* ones are in a planar geometry. The two next entries
* are undefined in that case.
*/
- bool isplanar = false;
- //!cnumber of connected hydrogens.
- int nHydrogens;
+ bool isplanar = false;
+ //! cnumber of connected hydrogens.
+ int nHydrogens;
//! Type for the heavy atom bonded to XH2/XH3.
std::string nextHeavyType;
//! The type of MNH* or MCH3* dummy mass to use.
*
* \param[in] name Residue name.
*/
- explicit VirtualSiteTopology(const std::string &name) : resname(name)
- {}
+ explicit VirtualSiteTopology(const std::string& name) : resname(name) {}
//! Residue name.
std::string resname;
//! Helper struct for single bond in virtual site.
* \param[in] a2 Second atom name.
* \param[in] v Value for distance.
*/
- VirtualSiteBond(const std::string &a1, const std::string &a2, real v) :
- atom1(a1), atom2(a2), value(v)
- {}
+ VirtualSiteBond(const std::string& a1, const std::string& a2, real v) :
+ atom1(a1),
+ atom2(a2),
+ value(v)
+ {
+ }
//! Atom 1 in bond.
std::string atom1;
//! Atom 2 in bond.
std::string atom2;
//! Distance value between atoms.
- float value;
+ float value;
};
//! Container of all bonds in virtual site.
std::vector<VirtualSiteBond> bond;
* \param[in] a3 Third atom name.
* \param[in] v Value for angle.
*/
- VirtualSiteAngle(const std::string &a1, const std::string &a2, const std::string &a3, real v) :
- atom1(a1), atom2(a2), atom3(a3), value(v)
- {}
+ VirtualSiteAngle(const std::string& a1, const std::string& a2, const std::string& a3, real v) :
+ atom1(a1),
+ atom2(a2),
+ atom3(a3),
+ value(v)
+ {
+ }
//! Atom 1 in angle.
std::string atom1;
//! Atom 2 in angle.
//! Atom 3 in angle.
std::string atom3;
//! Value for angle.
- float value;
+ float value;
};
//! Container for all angles in virtual site.
std::vector<VirtualSiteAngle> angle;
};
-enum {
- DDB_CH3, DDB_NH3, DDB_NH2, DDB_PHE, DDB_TYR,
- DDB_TRP, DDB_HISA, DDB_HISB, DDB_HISH, DDB_DIR_NR
+enum
+{
+ DDB_CH3,
+ DDB_NH3,
+ DDB_NH2,
+ DDB_PHE,
+ DDB_TYR,
+ DDB_TRP,
+ DDB_HISA,
+ DDB_HISB,
+ DDB_HISH,
+ DDB_DIR_NR
};
typedef char t_dirname[STRLEN];
-static const t_dirname ddb_dirnames[DDB_DIR_NR] = {
- "CH3",
- "NH3",
- "NH2",
- "PHE",
- "TYR",
- "TRP",
- "HISA",
- "HISB",
- "HISH"
-};
+static const t_dirname ddb_dirnames[DDB_DIR_NR] = { "CH3", "NH3", "NH2", "PHE", "TYR",
+ "TRP", "HISA", "HISB", "HISH" };
-static int ddb_name2dir(char *name)
+static int ddb_name2dir(char* name)
{
/* Translate a directive name to the number of the directive.
* HID/HIE/HIP names are translated to the ones we use in Gromacs.
}
-static void read_vsite_database(const char *ddbname,
- std::vector<VirtualSiteConfiguration> *vsiteconflist,
- std::vector<VirtualSiteTopology> *vsitetoplist)
+static void read_vsite_database(const char* ddbname,
+ std::vector<VirtualSiteConfiguration>* vsiteconflist,
+ std::vector<VirtualSiteTopology>* vsitetoplist)
{
/* This routine is a quick hack to fix the problem with hardcoded atomtypes
* and aromatic vsite parameters by reading them from a ff???.vsd file.
* case the second field should just be the word planar.
*/
- char dirstr[STRLEN];
- char pline[STRLEN];
- int curdir;
- char *ch;
- char s1[MAXNAME], s2[MAXNAME], s3[MAXNAME], s4[MAXNAME];
+ char dirstr[STRLEN];
+ char pline[STRLEN];
+ int curdir;
+ char* ch;
+ char s1[MAXNAME], s2[MAXNAME], s3[MAXNAME], s4[MAXNAME];
gmx::FilePtr ddb = gmx::openLibraryFile(ddbname);
curdir = -1;
- while (fgets2(pline, STRLEN-2, ddb.get()) != nullptr)
+ while (fgets2(pline, STRLEN - 2, ddb.get()) != nullptr)
{
strip_comment(pline);
trim(pline);
{
if (pline[0] == OPENDIR)
{
- strncpy(dirstr, pline+1, STRLEN-2);
- if ((ch = strchr (dirstr, CLOSEDIR)) != nullptr)
+ strncpy(dirstr, pline + 1, STRLEN - 2);
+ if ((ch = strchr(dirstr, CLOSEDIR)) != nullptr)
{
(*ch) = 0;
}
- trim (dirstr);
+ trim(dirstr);
- if (!gmx_strcasecmp(dirstr, "HID") ||
- !gmx_strcasecmp(dirstr, "HISD"))
+ if (!gmx_strcasecmp(dirstr, "HID") || !gmx_strcasecmp(dirstr, "HISD"))
{
sprintf(dirstr, "HISA");
}
- else if (!gmx_strcasecmp(dirstr, "HIE") ||
- !gmx_strcasecmp(dirstr, "HISE"))
+ else if (!gmx_strcasecmp(dirstr, "HIE") || !gmx_strcasecmp(dirstr, "HISE"))
{
sprintf(dirstr, "HISB");
}
curdir = ddb_name2dir(dirstr);
if (curdir < 0)
{
- gmx_fatal(FARGS, "Invalid directive %s in vsite database %s",
- dirstr, ddbname);
+ gmx_fatal(FARGS, "Invalid directive %s in vsite database %s", dirstr, ddbname);
}
}
else
case DDB_HISB:
case DDB_HISH:
{
- const auto found = std::find_if(vsitetoplist->begin(), vsitetoplist->end(),
- [&dirstr](const auto &entry)
- { return gmx::equalCaseInsensitive(dirstr, entry.resname); });
+ const auto found = std::find_if(
+ vsitetoplist->begin(), vsitetoplist->end(), [&dirstr](const auto& entry) {
+ return gmx::equalCaseInsensitive(dirstr, entry.resname);
+ });
/* Allocate a new topology entry if this is a new residue */
if (found == vsitetoplist->end())
{
else if (numberOfSites == 4)
{
/* angle */
- vsitetoplist->back().angle.emplace_back(s1String, s2String, s3String, strtod(s4, nullptr));
+ vsitetoplist->back().angle.emplace_back(s1String, s2String, s3String,
+ strtod(s4, nullptr));
/* angle */
}
else
{
- gmx_fatal(FARGS, "Need 3 or 4 values to specify bond/angle values in %s: %s\n", ddbname, pline);
+ gmx_fatal(FARGS,
+ "Need 3 or 4 values to specify bond/angle values in %s: %s\n",
+ ddbname, pline);
}
}
break;
default:
- gmx_fatal(FARGS, "Didnt find a case for directive %s in read_vsite_database\n", dirstr);
+ gmx_fatal(FARGS,
+ "Didnt find a case for directive %s in read_vsite_database\n", dirstr);
}
}
}
}
static int nitrogen_is_planar(gmx::ArrayRef<const VirtualSiteConfiguration> vsiteconflist,
- const std::string &atomtype)
+ const std::string& atomtype)
{
/* Return 1 if atomtype exists in database list and is planar, 0 if not,
* and -1 if not found.
*/
int res;
- const auto found = std::find_if(vsiteconflist.begin(), vsiteconflist.end(),
- [&atomtype](const auto &entry)
- { return (gmx::equalCaseInsensitive(entry.atomtype, atomtype) &&
- entry.nHydrogens == 2); });
+ const auto found =
+ std::find_if(vsiteconflist.begin(), vsiteconflist.end(), [&atomtype](const auto& entry) {
+ return (gmx::equalCaseInsensitive(entry.atomtype, atomtype) && entry.nHydrogens == 2);
+ });
if (found != vsiteconflist.end())
{
res = static_cast<int>(found->isplanar);
}
static std::string get_dummymass_name(gmx::ArrayRef<const VirtualSiteConfiguration> vsiteconflist,
- const std::string &atom, const std::string &nextheavy)
+ const std::string& atom,
+ const std::string& nextheavy)
{
/* Return the dummy mass name if found, or NULL if not set in ddb database */
- const auto found = std::find_if(vsiteconflist.begin(), vsiteconflist.end(),
- [&atom, &nextheavy](const auto &entry)
- { return (gmx::equalCaseInsensitive(atom, entry.atomtype) &&
- gmx::equalCaseInsensitive(nextheavy, entry.nextHeavyType)); });
+ const auto found = std::find_if(
+ vsiteconflist.begin(), vsiteconflist.end(), [&atom, &nextheavy](const auto& entry) {
+ return (gmx::equalCaseInsensitive(atom, entry.atomtype)
+ && gmx::equalCaseInsensitive(nextheavy, entry.nextHeavyType));
+ });
if (found != vsiteconflist.end())
{
return found->dummyMass;
}
-
static real get_ddb_bond(gmx::ArrayRef<const VirtualSiteTopology> vsitetop,
- const std::string &res,
- const std::string &atom1,
- const std::string &atom2)
+ const std::string& res,
+ const std::string& atom1,
+ const std::string& atom2)
{
- const auto found = std::find_if(vsitetop.begin(), vsitetop.end(),
- [&res](const auto &entry)
- { return gmx::equalCaseInsensitive(res, entry.resname); });
+ const auto found = std::find_if(vsitetop.begin(), vsitetop.end(), [&res](const auto& entry) {
+ return gmx::equalCaseInsensitive(res, entry.resname);
+ });
if (found == vsitetop.end())
{
gmx_fatal(FARGS, "No vsite information for residue %s found in vsite database.\n", res.c_str());
}
- const auto foundBond = std::find_if(found->bond.begin(), found->bond.end(),
- [&atom1, &atom2](const auto &entry)
- { return ((atom1 == entry.atom1 && atom2 == entry.atom2) ||
- (atom1 == entry.atom2 && atom2 == entry.atom1)); });
+ const auto foundBond =
+ std::find_if(found->bond.begin(), found->bond.end(), [&atom1, &atom2](const auto& entry) {
+ return ((atom1 == entry.atom1 && atom2 == entry.atom2)
+ || (atom1 == entry.atom2 && atom2 == entry.atom1));
+ });
if (foundBond == found->bond.end())
{
- gmx_fatal(FARGS, "Couldnt find bond %s-%s for residue %s in vsite database.\n", atom1.c_str(), atom2.c_str(), res.c_str());
+ gmx_fatal(FARGS, "Couldnt find bond %s-%s for residue %s in vsite database.\n",
+ atom1.c_str(), atom2.c_str(), res.c_str());
}
return foundBond->value;
static real get_ddb_angle(gmx::ArrayRef<const VirtualSiteTopology> vsitetop,
- const std::string &res,
- const std::string &atom1,
- const std::string &atom2,
- const std::string &atom3)
+ const std::string& res,
+ const std::string& atom1,
+ const std::string& atom2,
+ const std::string& atom3)
{
- const auto found = std::find_if(vsitetop.begin(), vsitetop.end(),
- [&res](const auto &entry)
- { return gmx::equalCaseInsensitive(res, entry.resname); });
+ const auto found = std::find_if(vsitetop.begin(), vsitetop.end(), [&res](const auto& entry) {
+ return gmx::equalCaseInsensitive(res, entry.resname);
+ });
if (found == vsitetop.end())
{
gmx_fatal(FARGS, "No vsite information for residue %s found in vsite database.\n", res.c_str());
}
- const auto foundAngle = std::find_if(found->angle.begin(), found->angle.end(),
- [&atom1, &atom2, &atom3](const auto &entry)
- { return ((atom1 == entry.atom1 && atom2 == entry.atom2 && atom3 == entry.atom3) ||
- (atom1 == entry.atom3 && atom2 == entry.atom2 && atom3 == entry.atom1) ||
- (atom1 == entry.atom2 && atom2 == entry.atom1 && atom3 == entry.atom3) ||
- (atom1 == entry.atom3 && atom2 == entry.atom1 && atom3 == entry.atom2)); });
+ const auto foundAngle = std::find_if(
+ found->angle.begin(), found->angle.end(), [&atom1, &atom2, &atom3](const auto& entry) {
+ return ((atom1 == entry.atom1 && atom2 == entry.atom2 && atom3 == entry.atom3)
+ || (atom1 == entry.atom3 && atom2 == entry.atom2 && atom3 == entry.atom1)
+ || (atom1 == entry.atom2 && atom2 == entry.atom1 && atom3 == entry.atom3)
+ || (atom1 == entry.atom3 && atom2 == entry.atom1 && atom3 == entry.atom2));
+ });
if (foundAngle == found->angle.end())
{
- 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());
+ 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());
}
return foundAngle->value;
}
-static void count_bonds(int atom, InteractionsOfType *psb, char ***atomname,
- int *nrbonds, int *nrHatoms, int Hatoms[], int *Heavy,
- int *nrheavies, int heavies[])
+static void count_bonds(int atom,
+ InteractionsOfType* psb,
+ char*** atomname,
+ int* nrbonds,
+ int* nrHatoms,
+ int Hatoms[],
+ int* Heavy,
+ int* nrheavies,
+ int heavies[])
{
int heavy, other, nrb, nrH, nrhv;
/* find heavy atom bound to this hydrogen */
heavy = NOTSET;
- for (auto parm = psb->interactionTypes.begin(); (parm != psb->interactionTypes.end()) && (heavy == NOTSET); parm++)
+ for (auto parm = psb->interactionTypes.begin();
+ (parm != psb->interactionTypes.end()) && (heavy == NOTSET); parm++)
{
if (parm->ai() == atom)
{
}
if (heavy == NOTSET)
{
- gmx_fatal(FARGS, "unbound hydrogen atom %d", atom+1);
+ gmx_fatal(FARGS, "unbound hydrogen atom %d", atom + 1);
}
/* find all atoms bound to heavy atom */
other = NOTSET;
nrb = 0;
nrH = 0;
nrhv = 0;
- for (const auto &parm : psb->interactionTypes)
+ for (const auto& parm : psb->interactionTypes)
{
if (parm.ai() == heavy)
{
*nrheavies = nrhv;
}
-static void print_bonds(FILE *fp, int o2n[],
- int nrHatoms, const int Hatoms[], int Heavy,
- int nrheavies, const int heavies[])
+static void
+print_bonds(FILE* fp, int o2n[], int nrHatoms, const int Hatoms[], int Heavy, int nrheavies, const int heavies[])
{
int i;
fprintf(fp, "Found: %d Hatoms: ", nrHatoms);
for (i = 0; i < nrHatoms; i++)
{
- fprintf(fp, " %d", o2n[Hatoms[i]]+1);
+ fprintf(fp, " %d", o2n[Hatoms[i]] + 1);
}
- fprintf(fp, "; %d Heavy atoms: %d", nrheavies+1, o2n[Heavy]+1);
+ fprintf(fp, "; %d Heavy atoms: %d", nrheavies + 1, o2n[Heavy] + 1);
for (i = 0; i < nrheavies; i++)
{
- fprintf(fp, " %d", o2n[heavies[i]]+1);
+ fprintf(fp, " %d", o2n[heavies[i]] + 1);
}
fprintf(fp, "\n");
}
-static int get_atype(int atom, t_atoms *at, gmx::ArrayRef<const PreprocessResidue> rtpFFDB,
- ResidueType *rt)
+static int get_atype(int atom, t_atoms* at, gmx::ArrayRef<const PreprocessResidue> rtpFFDB, ResidueType* rt)
{
- int type;
- bool bNterm;
+ int type;
+ bool bNterm;
if (at->atom[atom].m != 0.0F)
{
else
{
/* get type from rtpFFDB */
- auto localPpResidue = getDatabaseEntry(*(at->resinfo[at->atom[atom].resind].name), rtpFFDB);
- bNterm = rt->namedResidueHasType(*(at->resinfo[at->atom[atom].resind].name), "Protein") &&
- (at->atom[atom].resind == 0);
- int j = search_jtype(*localPpResidue, *(at->atomname[atom]), bNterm);
- type = localPpResidue->atom[j].type;
+ auto localPpResidue = getDatabaseEntry(*(at->resinfo[at->atom[atom].resind].name), rtpFFDB);
+ bNterm = rt->namedResidueHasType(*(at->resinfo[at->atom[atom].resind].name), "Protein")
+ && (at->atom[atom].resind == 0);
+ int j = search_jtype(*localPpResidue, *(at->atomname[atom]), bNterm);
+ type = localPpResidue->atom[j].type;
}
return type;
}
-static int vsite_nm2type(const char *name, PreprocessingAtomTypes *atype)
+static int vsite_nm2type(const char* name, PreprocessingAtomTypes* atype)
{
int tp;
tp = atype->atomTypeFromName(name);
if (tp == NOTSET)
{
- gmx_fatal(FARGS, "Dummy mass type (%s) not found in atom type database",
- name);
+ gmx_fatal(FARGS, "Dummy mass type (%s) not found in atom type database", name);
}
return tp;
}
-static real get_amass(int atom, t_atoms *at, gmx::ArrayRef<const PreprocessResidue> rtpFFDB,
- ResidueType *rt)
+static real get_amass(int atom, t_atoms* at, gmx::ArrayRef<const PreprocessResidue> rtpFFDB, ResidueType* rt)
{
- real mass;
- bool bNterm;
+ real mass;
+ bool bNterm;
if (at->atom[atom].m != 0.0F)
{
else
{
/* get mass from rtpFFDB */
- auto localPpResidue = getDatabaseEntry(*(at->resinfo[at->atom[atom].resind].name), rtpFFDB);
- bNterm = rt->namedResidueHasType(*(at->resinfo[at->atom[atom].resind].name), "Protein") &&
- (at->atom[atom].resind == 0);
- int j = search_jtype(*localPpResidue, *(at->atomname[atom]), bNterm);
- mass = localPpResidue->atom[j].m;
+ auto localPpResidue = getDatabaseEntry(*(at->resinfo[at->atom[atom].resind].name), rtpFFDB);
+ bNterm = rt->namedResidueHasType(*(at->resinfo[at->atom[atom].resind].name), "Protein")
+ && (at->atom[atom].resind == 0);
+ int j = search_jtype(*localPpResidue, *(at->atomname[atom]), bNterm);
+ mass = localPpResidue->atom[j].m;
}
return mass;
}
-static void my_add_param(InteractionsOfType *plist, int ai, int aj, real b)
+static void my_add_param(InteractionsOfType* plist, int ai, int aj, real b)
{
- static real c[MAXFORCEPARAM] =
- { NOTSET, NOTSET, NOTSET, NOTSET, NOTSET, NOTSET };
+ static real c[MAXFORCEPARAM] = { NOTSET, NOTSET, NOTSET, NOTSET, NOTSET, NOTSET };
c[0] = b;
add_param(plist, ai, aj, c, nullptr);
}
-static void add_vsites(gmx::ArrayRef<InteractionsOfType> plist, int vsite_type[],
- int Heavy, int nrHatoms, int Hatoms[],
- int nrheavies, int heavies[])
+static void add_vsites(gmx::ArrayRef<InteractionsOfType> plist,
+ int vsite_type[],
+ int Heavy,
+ int nrHatoms,
+ int Hatoms[],
+ int nrheavies,
+ int heavies[])
{
- int other, moreheavy;
+ int other, moreheavy;
for (int i = 0; i < nrHatoms; i++)
{
{
gmx_incons("Undetected error in setting up virtual sites");
}
- bool bSwapParity = (ftype < 0);
+ bool bSwapParity = (ftype < 0);
vsite_type[Hatoms[i]] = ftype = abs(ftype);
if (ftype == F_BONDS)
{
- if ( (nrheavies != 1) && (nrHatoms != 1) )
+ if ((nrheavies != 1) && (nrHatoms != 1))
{
- gmx_fatal(FARGS, "cannot make constraint in add_vsites for %d heavy "
- "atoms and %d hydrogen atoms", nrheavies, nrHatoms);
+ gmx_fatal(FARGS,
+ "cannot make constraint in add_vsites for %d heavy "
+ "atoms and %d hydrogen atoms",
+ nrheavies, nrHatoms);
}
my_add_param(&(plist[F_CONSTRNC]), Hatoms[i], heavies[0], NOTSET);
}
if (nrheavies < 2)
{
gmx_fatal(FARGS, "Not enough heavy atoms (%d) for %s (min 3)",
- nrheavies+1,
- interaction_function[vsite_type[Hatoms[i]]].name);
+ nrheavies + 1, interaction_function[vsite_type[Hatoms[i]]].name);
}
- add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], heavies[1],
- bSwapParity);
+ add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], heavies[1], bSwapParity);
break;
case F_VSITE3FAD:
{
/* find more heavy atoms */
other = moreheavy = NOTSET;
for (auto parm = plist[F_BONDS].interactionTypes.begin();
- (parm != plist[F_BONDS].interactionTypes.end()) && (moreheavy == NOTSET); parm++)
+ (parm != plist[F_BONDS].interactionTypes.end()) && (moreheavy == NOTSET);
+ parm++)
{
if (parm->ai() == heavies[0])
{
{
other = parm->ai();
}
- if ( (other != NOTSET) && (other != Heavy) )
+ if ((other != NOTSET) && (other != Heavy))
{
moreheavy = other;
}
}
if (moreheavy == NOTSET)
{
- gmx_fatal(FARGS, "Unbound molecule part %d-%d", Heavy+1, Hatoms[0]+1);
+ gmx_fatal(FARGS, "Unbound molecule part %d-%d", Heavy + 1, Hatoms[0] + 1);
}
}
- add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], moreheavy,
- bSwapParity);
+ add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], moreheavy, bSwapParity);
break;
}
case F_VSITE4FD:
if (nrheavies < 3)
{
gmx_fatal(FARGS, "Not enough heavy atoms (%d) for %s (min 4)",
- nrheavies+1,
- interaction_function[vsite_type[Hatoms[i]]].name);
+ nrheavies + 1, interaction_function[vsite_type[Hatoms[i]]].name);
}
- add_vsite4_atoms(&plist[ftype],
- Hatoms[0], Heavy, heavies[0], heavies[1], heavies[2]);
+ add_vsite4_atoms(&plist[ftype], Hatoms[0], Heavy, heavies[0], heavies[1], heavies[2]);
break;
default:
} /* for i */
}
-#define ANGLE_6RING (DEG2RAD*120)
+#define ANGLE_6RING (DEG2RAD * 120)
/* cosine rule: a^2 = b^2 + c^2 - 2 b c cos(alpha) */
/* get a^2 when a, b and alpha are given: */
-#define cosrule(b, c, alpha) ( gmx::square(b) + gmx::square(c) - 2*(b)*(c)*std::cos(alpha) )
+#define cosrule(b, c, alpha) (gmx::square(b) + gmx::square(c) - 2 * (b) * (c)*std::cos(alpha))
/* get cos(alpha) when a, b and c are given: */
-#define acosrule(a, b, c) ( (gmx::square(b)+gmx::square(c)-gmx::square(a))/(2*(b)*(c)) )
-
-static int gen_vsites_6ring(t_atoms *at, int *vsite_type[], gmx::ArrayRef<InteractionsOfType> plist,
- int nrfound, int *ats, real bond_cc, real bond_ch,
- real xcom, bool bDoZ)
+#define acosrule(a, b, c) ((gmx::square(b) + gmx::square(c) - gmx::square(a)) / (2 * (b) * (c)))
+
+static int gen_vsites_6ring(t_atoms* at,
+ int* vsite_type[],
+ gmx::ArrayRef<InteractionsOfType> plist,
+ int nrfound,
+ int* ats,
+ real bond_cc,
+ real bond_ch,
+ real xcom,
+ bool bDoZ)
{
/* these MUST correspond to the atnms array in do_vsite_aromatics! */
- enum {
- atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
- atCZ, atHZ, atNR
+ enum
+ {
+ atCG,
+ atCD1,
+ atHD1,
+ atCD2,
+ atHD2,
+ atCE1,
+ atHE1,
+ atCE2,
+ atHE2,
+ atCZ,
+ atHZ,
+ atNR
};
int i, nvsite;
}
/* constraints between CG, CE1 and CE2: */
- dCGCE = std::sqrt( cosrule(bond_cc, bond_cc, ANGLE_6RING) );
+ dCGCE = std::sqrt(cosrule(bond_cc, bond_cc, ANGLE_6RING));
my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE);
my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE2], dCGCE);
my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atCE2], dCGCE);
/* rest will be vsite3 */
mtot = 0;
nvsite = 0;
- for (i = 0; i < (bDoZ ? atNR : atHZ); i++)
+ for (i = 0; i < (bDoZ ? atNR : atHZ); i++)
{
mtot += at->atom[ats[i]].m;
- if (i != atCG && i != atCE1 && i != atCE2 && (bDoZ || (i != atHZ && i != atCZ) ) )
+ if (i != atCG && i != atCE1 && i != atCE2 && (bDoZ || (i != atHZ && i != atCZ)))
{
- at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
- (*vsite_type)[ats[i]] = F_VSITE3;
+ at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
+ (*vsite_type)[ats[i]] = F_VSITE3;
nvsite++;
}
}
* The center-of-mass in the call is defined with x=0 at
* the CE1-CE2 bond and y=0 at the line from CG to the middle of CE1-CE2 bond.
*/
- xCG = -bond_cc+bond_cc*std::cos(ANGLE_6RING);
+ xCG = -bond_cc + bond_cc * std::cos(ANGLE_6RING);
- mG = at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = xcom*mtot/xCG;
- mrest = mtot-mG;
- at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB =
- at->atom[ats[atCE2]].m = at->atom[ats[atCE2]].mB = mrest / 2;
+ mG = at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = xcom * mtot / xCG;
+ mrest = mtot - mG;
+ at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB = at->atom[ats[atCE2]].m =
+ at->atom[ats[atCE2]].mB = mrest / 2;
/* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
- tmp1 = dCGCE*std::sin(ANGLE_6RING*0.5);
- tmp2 = bond_cc*std::cos(0.5*ANGLE_6RING) + tmp1;
+ tmp1 = dCGCE * std::sin(ANGLE_6RING * 0.5);
+ tmp2 = bond_cc * std::cos(0.5 * ANGLE_6RING) + tmp1;
tmp1 *= 2;
- a = b = -bond_ch / tmp1;
+ a = b = -bond_ch / tmp1;
/* HE1 and HE2: */
- add_vsite3_param(&plist[F_VSITE3],
- ats[atHE1], ats[atCE1], ats[atCE2], ats[atCG], a, b);
- add_vsite3_param(&plist[F_VSITE3],
- ats[atHE2], ats[atCE2], ats[atCE1], ats[atCG], a, b);
+ add_vsite3_param(&plist[F_VSITE3], ats[atHE1], ats[atCE1], ats[atCE2], ats[atCG], a, b);
+ add_vsite3_param(&plist[F_VSITE3], ats[atHE2], ats[atCE2], ats[atCE1], ats[atCG], a, b);
/* CD1, CD2 and CZ: */
a = b = tmp2 / tmp1;
- add_vsite3_param(&plist[F_VSITE3],
- ats[atCD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
- add_vsite3_param(&plist[F_VSITE3],
- ats[atCD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
+ add_vsite3_param(&plist[F_VSITE3], ats[atCD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
+ add_vsite3_param(&plist[F_VSITE3], ats[atCD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
if (bDoZ)
{
- add_vsite3_param(&plist[F_VSITE3],
- ats[atCZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
+ add_vsite3_param(&plist[F_VSITE3], ats[atCZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
}
/* HD1, HD2 and HZ: */
- a = b = ( bond_ch + tmp2 ) / tmp1;
- add_vsite3_param(&plist[F_VSITE3],
- ats[atHD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
- add_vsite3_param(&plist[F_VSITE3],
- ats[atHD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
+ a = b = (bond_ch + tmp2) / tmp1;
+ add_vsite3_param(&plist[F_VSITE3], ats[atHD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
+ add_vsite3_param(&plist[F_VSITE3], ats[atHD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
if (bDoZ)
{
- add_vsite3_param(&plist[F_VSITE3],
- ats[atHZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
+ add_vsite3_param(&plist[F_VSITE3], ats[atHZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
}
return nvsite;
}
-static int gen_vsites_phe(t_atoms *at, int *vsite_type[], gmx::ArrayRef<InteractionsOfType> plist,
- int nrfound, int *ats, gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
+static int gen_vsites_phe(t_atoms* at,
+ int* vsite_type[],
+ gmx::ArrayRef<InteractionsOfType> plist,
+ int nrfound,
+ int* ats,
+ gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
{
real bond_cc, bond_ch;
real xcom, mtot;
int i;
/* these MUST correspond to the atnms array in do_vsite_aromatics! */
- enum {
- atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
- atCZ, atHZ, atNR
+ enum
+ {
+ atCG,
+ atCD1,
+ atHD1,
+ atCD2,
+ atHD2,
+ atCE1,
+ atHE1,
+ atCE2,
+ atHE2,
+ atCZ,
+ atHZ,
+ atNR
};
real x[atNR];
/* Aromatic rings have 6-fold symmetry, so we only need one bond length.
bond_cc = get_ddb_bond(vsitetop, "PHE", "CD1", "CE1");
bond_ch = get_ddb_bond(vsitetop, "PHE", "CD1", "HD1");
- x[atCG] = -bond_cc+bond_cc*std::cos(ANGLE_6RING);
+ x[atCG] = -bond_cc + bond_cc * std::cos(ANGLE_6RING);
x[atCD1] = -bond_cc;
- x[atHD1] = x[atCD1]+bond_ch*std::cos(ANGLE_6RING);
+ x[atHD1] = x[atCD1] + bond_ch * std::cos(ANGLE_6RING);
x[atCE1] = 0;
- x[atHE1] = x[atCE1]-bond_ch*std::cos(ANGLE_6RING);
+ x[atHE1] = x[atCE1] - bond_ch * std::cos(ANGLE_6RING);
x[atCD2] = x[atCD1];
x[atHD2] = x[atHD1];
x[atCE2] = x[atCE1];
x[atHE2] = x[atHE1];
- x[atCZ] = bond_cc*std::cos(0.5*ANGLE_6RING);
- x[atHZ] = x[atCZ]+bond_ch;
+ x[atCZ] = bond_cc * std::cos(0.5 * ANGLE_6RING);
+ x[atHZ] = x[atCZ] + bond_ch;
xcom = mtot = 0;
for (i = 0; i < atNR; i++)
{
- xcom += x[i]*at->atom[ats[i]].m;
+ xcom += x[i] * at->atom[ats[i]].m;
mtot += at->atom[ats[i]].m;
}
xcom /= mtot;
return gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, TRUE);
}
-static void calc_vsite3_param(real xd, real yd, real xi, real yi, real xj, real yj,
- real xk, real yk, real *a, real *b)
+static void
+calc_vsite3_param(real xd, real yd, real xi, real yi, real xj, real yj, real xk, real yk, real* a, real* b)
{
/* determine parameters by solving the equation system, since we know the
* virtual site coordinates here.
*/
real dx_ij, dx_ik, dy_ij, dy_ik;
- dx_ij = xj-xi;
- dy_ij = yj-yi;
- dx_ik = xk-xi;
- dy_ik = yk-yi;
+ dx_ij = xj - xi;
+ dy_ij = yj - yi;
+ dx_ik = xk - xi;
+ dy_ik = yk - yi;
- *a = ( (xd-xi)*dy_ik - dx_ik*(yd-yi) ) / (dx_ij*dy_ik - dx_ik*dy_ij);
- *b = ( yd - yi - (*a)*dy_ij ) / dy_ik;
+ *a = ((xd - xi) * dy_ik - dx_ik * (yd - yi)) / (dx_ij * dy_ik - dx_ik * dy_ij);
+ *b = (yd - yi - (*a) * dy_ij) / dy_ik;
}
-static int gen_vsites_trp(PreprocessingAtomTypes *atype,
- std::vector<gmx::RVec> *newx,
- t_atom *newatom[], char ***newatomname[],
- int *o2n[], int *newvsite_type[], int *newcgnr[],
- t_symtab *symtab, int *nadd,
- gmx::ArrayRef<const gmx::RVec> x, int *cgnr[],
- t_atoms *at, int *vsite_type[],
- gmx::ArrayRef<InteractionsOfType> plist,
- int nrfound, int *ats, int add_shift,
+static int gen_vsites_trp(PreprocessingAtomTypes* atype,
+ std::vector<gmx::RVec>* newx,
+ t_atom* newatom[],
+ char*** newatomname[],
+ int* o2n[],
+ int* newvsite_type[],
+ int* newcgnr[],
+ t_symtab* symtab,
+ int* nadd,
+ gmx::ArrayRef<const gmx::RVec> x,
+ int* cgnr[],
+ t_atoms* at,
+ int* vsite_type[],
+ gmx::ArrayRef<InteractionsOfType> plist,
+ int nrfound,
+ int* ats,
+ int add_shift,
gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
{
#define NMASS 2
/* these MUST correspond to the atnms array in do_vsite_aromatics! */
- enum {
- atCB, atCG, atCD1, atHD1, atCD2, atNE1, atHE1, atCE2, atCE3, atHE3,
- atCZ2, atHZ2, atCZ3, atHZ3, atCH2, atHH2, atNR
+ enum
+ {
+ atCB,
+ atCG,
+ atCD1,
+ atHD1,
+ atCD2,
+ atNE1,
+ atHE1,
+ atCE2,
+ atCE3,
+ atHE3,
+ atCZ2,
+ atHZ2,
+ atCZ3,
+ atHZ3,
+ atCH2,
+ atHH2,
+ atNR
};
/* weights for determining the COM's of both rings (M1 and M2): */
- real mw[NMASS][atNR] = {
- { 0, 1, 1, 1, 0.5, 1, 1, 0.5, 0, 0,
- 0, 0, 0, 0, 0, 0 },
- { 0, 0, 0, 0, 0.5, 0, 0, 0.5, 1, 1,
- 1, 1, 1, 1, 1, 1 }
- };
+ real mw[NMASS][atNR] = { { 0, 1, 1, 1, 0.5, 1, 1, 0.5, 0, 0, 0, 0, 0, 0, 0, 0 },
+ { 0, 0, 0, 0, 0.5, 0, 0, 0.5, 1, 1, 1, 1, 1, 1, 1, 1 } };
real xi[atNR], yi[atNR];
real xcom[NMASS], ycom[NMASS], alpha;
b_CE3_HE3 = get_ddb_bond(vsitetop, "TRP", "CE3", "HE3");
b_CZ3_HZ3 = get_ddb_bond(vsitetop, "TRP", "CZ3", "HZ3");
- a_NE1_CE2_CD2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "NE1", "CE2", "CD2");
- a_CE2_CD2_CG = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CE2", "CD2", "CG");
- a_CB_CG_CD2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CB", "CG", "CD2");
- a_CD2_CG_CD1 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CD2", "CG", "CD1");
-
- a_CE2_CD2_CE3 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CE2", "CD2", "CE3");
- a_CD2_CE2_CZ2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CD2", "CE2", "CZ2");
- a_CD2_CE3_CZ3 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CD2", "CE3", "CZ3");
- a_CE3_CZ3_HZ3 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CE3", "CZ3", "HZ3");
- a_CZ2_CH2_HH2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CZ2", "CH2", "HH2");
- a_CE2_CZ2_HZ2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CE2", "CZ2", "HZ2");
- a_CE2_CZ2_CH2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CE2", "CZ2", "CH2");
- a_CG_CD1_HD1 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CG", "CD1", "HD1");
- a_HE1_NE1_CE2 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "HE1", "NE1", "CE2");
- a_CD2_CE3_HE3 = DEG2RAD*get_ddb_angle(vsitetop, "TRP", "CD2", "CE3", "HE3");
+ a_NE1_CE2_CD2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "NE1", "CE2", "CD2");
+ a_CE2_CD2_CG = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE2", "CD2", "CG");
+ a_CB_CG_CD2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CB", "CG", "CD2");
+ a_CD2_CG_CD1 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CD2", "CG", "CD1");
+
+ a_CE2_CD2_CE3 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE2", "CD2", "CE3");
+ a_CD2_CE2_CZ2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CD2", "CE2", "CZ2");
+ a_CD2_CE3_CZ3 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CD2", "CE3", "CZ3");
+ a_CE3_CZ3_HZ3 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE3", "CZ3", "HZ3");
+ a_CZ2_CH2_HH2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CZ2", "CH2", "HH2");
+ a_CE2_CZ2_HZ2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE2", "CZ2", "HZ2");
+ a_CE2_CZ2_CH2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CE2", "CZ2", "CH2");
+ a_CG_CD1_HD1 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CG", "CD1", "HD1");
+ a_HE1_NE1_CE2 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "HE1", "NE1", "CE2");
+ a_CD2_CE3_HE3 = DEG2RAD * get_ddb_angle(vsitetop, "TRP", "CD2", "CE3", "HE3");
/* Calculate local coordinates.
* y-axis (x=0) is the bond CD2-CE2.
* intersects the middle of the bond.
*/
xi[atCD2] = 0;
- yi[atCD2] = -0.5*b_CD2_CE2;
+ yi[atCD2] = -0.5 * b_CD2_CE2;
xi[atCE2] = 0;
- yi[atCE2] = 0.5*b_CD2_CE2;
+ yi[atCE2] = 0.5 * b_CD2_CE2;
- xi[atNE1] = -b_NE1_CE2*std::sin(a_NE1_CE2_CD2);
- yi[atNE1] = yi[atCE2]-b_NE1_CE2*std::cos(a_NE1_CE2_CD2);
+ xi[atNE1] = -b_NE1_CE2 * std::sin(a_NE1_CE2_CD2);
+ yi[atNE1] = yi[atCE2] - b_NE1_CE2 * std::cos(a_NE1_CE2_CD2);
- xi[atCG] = -b_CG_CD2*std::sin(a_CE2_CD2_CG);
- yi[atCG] = yi[atCD2]+b_CG_CD2*std::cos(a_CE2_CD2_CG);
+ xi[atCG] = -b_CG_CD2 * std::sin(a_CE2_CD2_CG);
+ yi[atCG] = yi[atCD2] + b_CG_CD2 * std::cos(a_CE2_CD2_CG);
alpha = a_CE2_CD2_CG + M_PI - a_CB_CG_CD2;
- xi[atCB] = xi[atCG]-b_CB_CG*std::sin(alpha);
- yi[atCB] = yi[atCG]+b_CB_CG*std::cos(alpha);
+ xi[atCB] = xi[atCG] - b_CB_CG * std::sin(alpha);
+ yi[atCB] = yi[atCG] + b_CB_CG * std::cos(alpha);
alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - M_PI;
- xi[atCD1] = xi[atCG]-b_CG_CD1*std::sin(alpha);
- yi[atCD1] = yi[atCG]+b_CG_CD1*std::cos(alpha);
+ xi[atCD1] = xi[atCG] - b_CG_CD1 * std::sin(alpha);
+ yi[atCD1] = yi[atCG] + b_CG_CD1 * std::cos(alpha);
- xi[atCE3] = b_CD2_CE3*std::sin(a_CE2_CD2_CE3);
- yi[atCE3] = yi[atCD2]+b_CD2_CE3*std::cos(a_CE2_CD2_CE3);
+ xi[atCE3] = b_CD2_CE3 * std::sin(a_CE2_CD2_CE3);
+ yi[atCE3] = yi[atCD2] + b_CD2_CE3 * std::cos(a_CE2_CD2_CE3);
- xi[atCZ2] = b_CE2_CZ2*std::sin(a_CD2_CE2_CZ2);
- yi[atCZ2] = yi[atCE2]-b_CE2_CZ2*std::cos(a_CD2_CE2_CZ2);
+ xi[atCZ2] = b_CE2_CZ2 * std::sin(a_CD2_CE2_CZ2);
+ yi[atCZ2] = yi[atCE2] - b_CE2_CZ2 * std::cos(a_CD2_CE2_CZ2);
alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - M_PI;
- xi[atCZ3] = xi[atCE3]+b_CE3_CZ3*std::sin(alpha);
- yi[atCZ3] = yi[atCE3]+b_CE3_CZ3*std::cos(alpha);
+ xi[atCZ3] = xi[atCE3] + b_CE3_CZ3 * std::sin(alpha);
+ yi[atCZ3] = yi[atCE3] + b_CE3_CZ3 * std::cos(alpha);
alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - M_PI;
- xi[atCH2] = xi[atCZ2]+b_CZ2_CH2*std::sin(alpha);
- yi[atCH2] = yi[atCZ2]-b_CZ2_CH2*std::cos(alpha);
+ xi[atCH2] = xi[atCZ2] + b_CZ2_CH2 * std::sin(alpha);
+ yi[atCH2] = yi[atCZ2] - b_CZ2_CH2 * std::cos(alpha);
/* hydrogens */
alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - a_CG_CD1_HD1;
- xi[atHD1] = xi[atCD1]-b_CD1_HD1*std::sin(alpha);
- yi[atHD1] = yi[atCD1]+b_CD1_HD1*std::cos(alpha);
+ xi[atHD1] = xi[atCD1] - b_CD1_HD1 * std::sin(alpha);
+ yi[atHD1] = yi[atCD1] + b_CD1_HD1 * std::cos(alpha);
alpha = a_NE1_CE2_CD2 + M_PI - a_HE1_NE1_CE2;
- xi[atHE1] = xi[atNE1]-b_NE1_HE1*std::sin(alpha);
- yi[atHE1] = yi[atNE1]-b_NE1_HE1*std::cos(alpha);
+ xi[atHE1] = xi[atNE1] - b_NE1_HE1 * std::sin(alpha);
+ yi[atHE1] = yi[atNE1] - b_NE1_HE1 * std::cos(alpha);
alpha = a_CE2_CD2_CE3 + M_PI - a_CD2_CE3_HE3;
- xi[atHE3] = xi[atCE3]+b_CE3_HE3*std::sin(alpha);
- yi[atHE3] = yi[atCE3]+b_CE3_HE3*std::cos(alpha);
+ xi[atHE3] = xi[atCE3] + b_CE3_HE3 * std::sin(alpha);
+ yi[atHE3] = yi[atCE3] + b_CE3_HE3 * std::cos(alpha);
alpha = a_CD2_CE2_CZ2 + M_PI - a_CE2_CZ2_HZ2;
- xi[atHZ2] = xi[atCZ2]+b_CZ2_HZ2*std::sin(alpha);
- yi[atHZ2] = yi[atCZ2]-b_CZ2_HZ2*std::cos(alpha);
+ xi[atHZ2] = xi[atCZ2] + b_CZ2_HZ2 * std::sin(alpha);
+ yi[atHZ2] = yi[atCZ2] - b_CZ2_HZ2 * std::cos(alpha);
alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - a_CZ2_CH2_HH2;
- xi[atHZ3] = xi[atCZ3]+b_CZ3_HZ3*std::sin(alpha);
- yi[atHZ3] = yi[atCZ3]+b_CZ3_HZ3*std::cos(alpha);
+ xi[atHZ3] = xi[atCZ3] + b_CZ3_HZ3 * std::sin(alpha);
+ yi[atHZ3] = yi[atCZ3] + b_CZ3_HZ3 * std::cos(alpha);
alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - a_CE3_CZ3_HZ3;
- xi[atHH2] = xi[atCH2]+b_CH2_HH2*std::sin(alpha);
- yi[atHH2] = yi[atCH2]-b_CH2_HH2*std::cos(alpha);
+ xi[atHH2] = xi[atCH2] + b_CH2_HH2 * std::sin(alpha);
+ yi[atHH2] = yi[atCH2] - b_CH2_HH2 * std::cos(alpha);
/* Calculate masses for each ring and put it on the dummy masses */
for (j = 0; j < NMASS; j++)
{
for (j = 0; j < NMASS; j++)
{
- mM[j] += mw[j][i] * at->atom[ats[i]].m;
+ mM[j] += mw[j][i] * at->atom[ats[i]].m;
xcom[j] += xi[i] * mw[j][i] * at->atom[ats[i]].m;
ycom[j] += yi[i] * mw[j][i] * at->atom[ats[i]].m;
}
i0 = ats[atCB];
for (j = 0; j < NMASS; j++)
{
- atM[j] = i0+*nadd+j;
+ atM[j] = i0 + *nadd + j;
}
if (debug)
{
- fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, (*o2n)[i0]+1);
+ fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, (*o2n)[i0] + 1);
}
*nadd += NMASS;
for (j = i0; j < at->nr; j++)
{
- (*o2n)[j] = j+*nadd;
+ (*o2n)[j] = j + *nadd;
}
- newx->resize(at->nr+*nadd);
- srenew(*newatom, at->nr+*nadd);
- srenew(*newatomname, at->nr+*nadd);
- srenew(*newvsite_type, at->nr+*nadd);
- srenew(*newcgnr, at->nr+*nadd);
+ newx->resize(at->nr + *nadd);
+ srenew(*newatom, at->nr + *nadd);
+ srenew(*newatomname, at->nr + *nadd);
+ srenew(*newvsite_type, at->nr + *nadd);
+ srenew(*newcgnr, at->nr + *nadd);
for (j = 0; j < NMASS; j++)
{
- (*newatomname)[at->nr+*nadd-1-j] = nullptr;
+ (*newatomname)[at->nr + *nadd - 1 - j] = nullptr;
}
/* Dummy masses will be placed at the center-of-mass in each ring. */
*/
rvec_sub(x[ats[atCB]], x[ats[atCG]], r_ij);
rvec_sub(x[ats[atCD2]], x[ats[atCG]], r_ik);
- calc_vsite3_param(xcom[0], ycom[0], xi[atCG], yi[atCG], xi[atCB], yi[atCB],
- xi[atCD2], yi[atCD2], &a, &b);
+ calc_vsite3_param(xcom[0], ycom[0], xi[atCG], yi[atCG], xi[atCB], yi[atCB], xi[atCD2],
+ yi[atCD2], &a, &b);
svmul(a, r_ij, t1);
svmul(b, r_ik, t2);
rvec_add(t1, t2, t1);
rvec_add(t1, x[ats[atCG]], (*newx)[atM[0]]);
- calc_vsite3_param(xcom[1], ycom[1], xi[atCG], yi[atCG], xi[atCB], yi[atCB],
- xi[atCD2], yi[atCD2], &a, &b);
+ calc_vsite3_param(xcom[1], ycom[1], xi[atCG], yi[atCG], xi[atCB], yi[atCB], xi[atCD2],
+ yi[atCD2], &a, &b);
svmul(a, r_ij, t1);
svmul(b, r_ik, t2);
rvec_add(t1, t2, t1);
/* set parameters for the masses */
for (j = 0; j < NMASS; j++)
{
- sprintf(name, "MW%d", j+1);
- (*newatomname) [atM[j]] = put_symtab(symtab, name);
- (*newatom) [atM[j]].m = (*newatom)[atM[j]].mB = mM[j];
- (*newatom) [atM[j]].q = (*newatom)[atM[j]].qB = 0.0;
- (*newatom) [atM[j]].type = (*newatom)[atM[j]].typeB = tpM;
- (*newatom) [atM[j]].ptype = eptAtom;
- (*newatom) [atM[j]].resind = at->atom[i0].resind;
- (*newatom) [atM[j]].elem[0] = 'M';
- (*newatom) [atM[j]].elem[1] = '\0';
- (*newvsite_type)[atM[j]] = NOTSET;
- (*newcgnr) [atM[j]] = (*cgnr)[i0];
+ sprintf(name, "MW%d", j + 1);
+ (*newatomname)[atM[j]] = put_symtab(symtab, name);
+ (*newatom)[atM[j]].m = (*newatom)[atM[j]].mB = mM[j];
+ (*newatom)[atM[j]].q = (*newatom)[atM[j]].qB = 0.0;
+ (*newatom)[atM[j]].type = (*newatom)[atM[j]].typeB = tpM;
+ (*newatom)[atM[j]].ptype = eptAtom;
+ (*newatom)[atM[j]].resind = at->atom[i0].resind;
+ (*newatom)[atM[j]].elem[0] = 'M';
+ (*newatom)[atM[j]].elem[1] = '\0';
+ (*newvsite_type)[atM[j]] = NOTSET;
+ (*newcgnr)[atM[j]] = (*cgnr)[i0];
}
/* renumber cgnr: */
for (i = i0; i < at->nr; i++)
/* constraints between CB, M1 and M2 */
/* 'add_shift' says which atoms won't be renumbered afterwards */
- dCBM1 = std::hypot( xcom[0]-xi[atCB], ycom[0]-yi[atCB] );
- dM1M2 = std::hypot( xcom[0]-xcom[1], ycom[0]-ycom[1] );
- dCBM2 = std::hypot( xcom[1]-xi[atCB], ycom[1]-yi[atCB] );
- my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift+atM[0], dCBM1);
- my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift+atM[1], dCBM2);
- my_add_param(&(plist[F_CONSTRNC]), add_shift+atM[0], add_shift+atM[1], dM1M2);
+ dCBM1 = std::hypot(xcom[0] - xi[atCB], ycom[0] - yi[atCB]);
+ dM1M2 = std::hypot(xcom[0] - xcom[1], ycom[0] - ycom[1]);
+ dCBM2 = std::hypot(xcom[1] - xi[atCB], ycom[1] - yi[atCB]);
+ my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift + atM[0], dCBM1);
+ my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift + atM[1], dCBM2);
+ my_add_param(&(plist[F_CONSTRNC]), add_shift + atM[0], add_shift + atM[1], dM1M2);
/* rest will be vsite3 */
nvsite = 0;
{
if (i != atCB)
{
- at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
- (*vsite_type)[ats[i]] = F_VSITE3;
+ at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
+ (*vsite_type)[ats[i]] = F_VSITE3;
nvsite++;
}
}
r_d = r_M1 + a r_M1_M2 + b r_M1_CB */
for (i = 0; i < atNR; i++)
{
- if ( (*vsite_type)[ats[i]] == F_VSITE3)
+ if ((*vsite_type)[ats[i]] == F_VSITE3)
{
- calc_vsite3_param(xi[i], yi[i], xcom[0], ycom[0], xcom[1], ycom[1], xi[atCB], yi[atCB], &a, &b);
- add_vsite3_param(&plist[F_VSITE3],
- ats[i], add_shift+atM[0], add_shift+atM[1], ats[atCB], a, b);
+ calc_vsite3_param(xi[i], yi[i], xcom[0], ycom[0], xcom[1], ycom[1], xi[atCB], yi[atCB],
+ &a, &b);
+ add_vsite3_param(&plist[F_VSITE3], ats[i], add_shift + atM[0], add_shift + atM[1],
+ ats[atCB], a, b);
}
}
return nvsite;
}
-static int gen_vsites_tyr(PreprocessingAtomTypes *atype,
- std::vector<gmx::RVec> *newx,
- t_atom *newatom[], char ***newatomname[],
- int *o2n[], int *newvsite_type[], int *newcgnr[],
- t_symtab *symtab, int *nadd,
- gmx::ArrayRef<const gmx::RVec> x, int *cgnr[],
- t_atoms *at, int *vsite_type[],
- gmx::ArrayRef<InteractionsOfType> plist,
- int nrfound, int *ats, int add_shift,
+static int gen_vsites_tyr(PreprocessingAtomTypes* atype,
+ std::vector<gmx::RVec>* newx,
+ t_atom* newatom[],
+ char*** newatomname[],
+ int* o2n[],
+ int* newvsite_type[],
+ int* newcgnr[],
+ t_symtab* symtab,
+ int* nadd,
+ gmx::ArrayRef<const gmx::RVec> x,
+ int* cgnr[],
+ t_atoms* at,
+ int* vsite_type[],
+ gmx::ArrayRef<InteractionsOfType> plist,
+ int nrfound,
+ int* ats,
+ int add_shift,
gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
{
int nvsite, i, i0, j, atM, tpM;
char name[10];
/* these MUST correspond to the atnms array in do_vsite_aromatics! */
- enum {
- atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
- atCZ, atOH, atHH, atNR
+ enum
+ {
+ atCG,
+ atCD1,
+ atHD1,
+ atCD2,
+ atHD2,
+ atCE1,
+ atHE1,
+ atCE2,
+ atHE2,
+ atCZ,
+ atOH,
+ atHH,
+ atNR
};
real xi[atNR];
/* CG, CE1, CE2 (as in general 6-ring) and OH and HH stay,
bond_ch = get_ddb_bond(vsitetop, "TYR", "CD1", "HD1");
bond_co = get_ddb_bond(vsitetop, "TYR", "CZ", "OH");
bond_oh = get_ddb_bond(vsitetop, "TYR", "OH", "HH");
- angle_coh = DEG2RAD*get_ddb_angle(vsitetop, "TYR", "CZ", "OH", "HH");
+ angle_coh = DEG2RAD * get_ddb_angle(vsitetop, "TYR", "CZ", "OH", "HH");
- xi[atCG] = -bond_cc+bond_cc*std::cos(ANGLE_6RING);
+ xi[atCG] = -bond_cc + bond_cc * std::cos(ANGLE_6RING);
xi[atCD1] = -bond_cc;
- xi[atHD1] = xi[atCD1]+bond_ch*std::cos(ANGLE_6RING);
+ xi[atHD1] = xi[atCD1] + bond_ch * std::cos(ANGLE_6RING);
xi[atCE1] = 0;
- xi[atHE1] = xi[atCE1]-bond_ch*std::cos(ANGLE_6RING);
+ xi[atHE1] = xi[atCE1] - bond_ch * std::cos(ANGLE_6RING);
xi[atCD2] = xi[atCD1];
xi[atHD2] = xi[atHD1];
xi[atCE2] = xi[atCE1];
xi[atHE2] = xi[atHE1];
- xi[atCZ] = bond_cc*std::cos(0.5*ANGLE_6RING);
- xi[atOH] = xi[atCZ]+bond_co;
+ xi[atCZ] = bond_cc * std::cos(0.5 * ANGLE_6RING);
+ xi[atOH] = xi[atCZ] + bond_co;
xcom = mtot = 0;
for (i = 0; i < atOH; i++)
{
- xcom += xi[i]*at->atom[ats[i]].m;
+ xcom += xi[i] * at->atom[ats[i]].m;
mtot += at->atom[ats[i]].m;
}
xcom /= mtot;
/* then construct CZ from the 2nd triangle */
/* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
- a = b = 0.5 * bond_co / ( bond_co - bond_cc*std::cos(ANGLE_6RING) );
- add_vsite3_param(&plist[F_VSITE3],
- ats[atCZ], ats[atOH], ats[atCE1], ats[atCE2], a, b);
+ a = b = 0.5 * bond_co / (bond_co - bond_cc * std::cos(ANGLE_6RING));
+ add_vsite3_param(&plist[F_VSITE3], ats[atCZ], ats[atOH], ats[atCE1], ats[atCE2], a, b);
at->atom[ats[atCZ]].m = at->atom[ats[atCZ]].mB = 0;
/* constraints between CE1, CE2 and OH */
- dCGCE = std::sqrt( cosrule(bond_cc, bond_cc, ANGLE_6RING) );
- dCEOH = std::sqrt( cosrule(bond_cc, bond_co, ANGLE_6RING) );
+ dCGCE = std::sqrt(cosrule(bond_cc, bond_cc, ANGLE_6RING));
+ dCEOH = std::sqrt(cosrule(bond_cc, bond_co, ANGLE_6RING));
my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atOH], dCEOH);
my_add_param(&(plist[F_CONSTRNC]), ats[atCE2], ats[atOH], dCEOH);
* apply to the HH constructed atom and not directly on the virtual mass.
*/
- vdist = 2.0*bond_oh;
- mM = at->atom[ats[atHH]].m/2.0;
- at->atom[ats[atOH]].m += mM; /* add 1/2 of original H mass */
+ vdist = 2.0 * bond_oh;
+ mM = at->atom[ats[atHH]].m / 2.0;
+ at->atom[ats[atOH]].m += mM; /* add 1/2 of original H mass */
at->atom[ats[atOH]].mB += mM; /* add 1/2 of original H mass */
- at->atom[ats[atHH]].m = at->atom[ats[atHH]].mB = 0;
+ at->atom[ats[atHH]].m = at->atom[ats[atHH]].mB = 0;
/* get dummy mass type */
tpM = vsite_nm2type("MW", atype);
/* make space for 1 mass: shift HH only */
i0 = ats[atHH];
- atM = i0+*nadd;
+ atM = i0 + *nadd;
if (debug)
{
- fprintf(stderr, "Inserting 1 dummy mass at %d\n", (*o2n)[i0]+1);
+ fprintf(stderr, "Inserting 1 dummy mass at %d\n", (*o2n)[i0] + 1);
}
(*nadd)++;
for (j = i0; j < at->nr; j++)
{
- (*o2n)[j] = j+*nadd;
+ (*o2n)[j] = j + *nadd;
}
- newx->resize(at->nr+*nadd);
- srenew(*newatom, at->nr+*nadd);
- srenew(*newatomname, at->nr+*nadd);
- srenew(*newvsite_type, at->nr+*nadd);
- srenew(*newcgnr, at->nr+*nadd);
- (*newatomname)[at->nr+*nadd-1] = nullptr;
+ newx->resize(at->nr + *nadd);
+ srenew(*newatom, at->nr + *nadd);
+ srenew(*newatomname, at->nr + *nadd);
+ srenew(*newvsite_type, at->nr + *nadd);
+ srenew(*newcgnr, at->nr + *nadd);
+ (*newatomname)[at->nr + *nadd - 1] = nullptr;
/* Calc the dummy mass initial position */
rvec_sub(x[ats[atHH]], x[ats[atOH]], r1);
rvec_add(r1, x[ats[atHH]], (*newx)[atM]);
strcpy(name, "MW1");
- (*newatomname) [atM] = put_symtab(symtab, name);
- (*newatom) [atM].m = (*newatom)[atM].mB = mM;
- (*newatom) [atM].q = (*newatom)[atM].qB = 0.0;
- (*newatom) [atM].type = (*newatom)[atM].typeB = tpM;
- (*newatom) [atM].ptype = eptAtom;
- (*newatom) [atM].resind = at->atom[i0].resind;
- (*newatom) [atM].elem[0] = 'M';
- (*newatom) [atM].elem[1] = '\0';
- (*newvsite_type)[atM] = NOTSET;
- (*newcgnr) [atM] = (*cgnr)[i0];
+ (*newatomname)[atM] = put_symtab(symtab, name);
+ (*newatom)[atM].m = (*newatom)[atM].mB = mM;
+ (*newatom)[atM].q = (*newatom)[atM].qB = 0.0;
+ (*newatom)[atM].type = (*newatom)[atM].typeB = tpM;
+ (*newatom)[atM].ptype = eptAtom;
+ (*newatom)[atM].resind = at->atom[i0].resind;
+ (*newatom)[atM].elem[0] = 'M';
+ (*newatom)[atM].elem[1] = '\0';
+ (*newvsite_type)[atM] = NOTSET;
+ (*newcgnr)[atM] = (*cgnr)[i0];
/* renumber cgnr: */
for (i = i0; i < at->nr; i++)
{
(*vsite_type)[ats[atHH]] = F_VSITE2;
nvsite++;
/* assume we also want the COH angle constrained: */
- tmp1 = bond_cc*std::cos(0.5*ANGLE_6RING) + dCGCE*std::sin(ANGLE_6RING*0.5) + bond_co;
- dCGM = std::sqrt( cosrule(tmp1, vdist, angle_coh) );
- my_add_param(&(plist[F_CONSTRNC]), ats[atCG], add_shift+atM, dCGM);
- my_add_param(&(plist[F_CONSTRNC]), ats[atOH], add_shift+atM, vdist);
+ tmp1 = bond_cc * std::cos(0.5 * ANGLE_6RING) + dCGCE * std::sin(ANGLE_6RING * 0.5) + bond_co;
+ dCGM = std::sqrt(cosrule(tmp1, vdist, angle_coh));
+ my_add_param(&(plist[F_CONSTRNC]), ats[atCG], add_shift + atM, dCGM);
+ my_add_param(&(plist[F_CONSTRNC]), ats[atOH], add_shift + atM, vdist);
- add_vsite2_param(&plist[F_VSITE2],
- ats[atHH], ats[atOH], add_shift+atM, 1.0/2.0);
+ add_vsite2_param(&plist[F_VSITE2], ats[atHH], ats[atOH], add_shift + atM, 1.0 / 2.0);
return nvsite;
}
-static int gen_vsites_his(t_atoms *at, int *vsite_type[],
- gmx::ArrayRef<InteractionsOfType> plist,
- int nrfound, int *ats, gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
+static int gen_vsites_his(t_atoms* at,
+ int* vsite_type[],
+ gmx::ArrayRef<InteractionsOfType> plist,
+ int nrfound,
+ int* ats,
+ gmx::ArrayRef<const VirtualSiteTopology> vsitetop)
{
int nvsite, i;
real a, b, alpha, dCGCE1, dCGNE2;
char resname[10];
/* these MUST correspond to the atnms array in do_vsite_aromatics! */
- enum {
- atCG, atND1, atHD1, atCD2, atHD2, atCE1, atHE1, atNE2, atHE2, atNR
+ enum
+ {
+ atCG,
+ atND1,
+ atHD1,
+ atCD2,
+ atHD2,
+ atCE1,
+ atHE1,
+ atNE2,
+ atHE2,
+ atNR
};
real x[atNR], y[atNR];
/* assert( nrfound >= atNR-3 || nrfound <= atNR );
* Don't understand the above logic. Shouldn't it be && rather than || ???
*/
- if ((nrfound < atNR-3) || (nrfound > atNR))
+ if ((nrfound < atNR - 3) || (nrfound > atNR))
{
gmx_incons("Generating vsites for HIS");
}
/* avoid warnings about uninitialized variables */
- b_ND1_HD1 = b_NE2_HE2 = b_CE1_HE1 = b_CD2_HD2 = a_NE2_CE1_HE1 =
- a_NE2_CD2_HD2 = a_CE1_ND1_HD1 = a_CE1_NE2_HE2 = 0;
+ b_ND1_HD1 = b_NE2_HE2 = b_CE1_HE1 = b_CD2_HD2 = a_NE2_CE1_HE1 = a_NE2_CD2_HD2 = a_CE1_ND1_HD1 =
+ a_CE1_NE2_HE2 = 0;
if (ats[atHD1] != NOTSET)
{
b_CE1_NE2 = get_ddb_bond(vsitetop, resname, "CE1", "NE2");
b_CG_CD2 = get_ddb_bond(vsitetop, resname, "CG", "CD2");
b_CD2_NE2 = get_ddb_bond(vsitetop, resname, "CD2", "NE2");
- a_CG_ND1_CE1 = DEG2RAD*get_ddb_angle(vsitetop, resname, "CG", "ND1", "CE1");
- a_CG_CD2_NE2 = DEG2RAD*get_ddb_angle(vsitetop, resname, "CG", "CD2", "NE2");
- a_ND1_CE1_NE2 = DEG2RAD*get_ddb_angle(vsitetop, resname, "ND1", "CE1", "NE2");
- a_CE1_NE2_CD2 = DEG2RAD*get_ddb_angle(vsitetop, resname, "CE1", "NE2", "CD2");
+ a_CG_ND1_CE1 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CG", "ND1", "CE1");
+ a_CG_CD2_NE2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CG", "CD2", "NE2");
+ a_ND1_CE1_NE2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "ND1", "CE1", "NE2");
+ a_CE1_NE2_CD2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CE1", "NE2", "CD2");
if (ats[atHD1] != NOTSET)
{
b_ND1_HD1 = get_ddb_bond(vsitetop, resname, "ND1", "HD1");
- a_CE1_ND1_HD1 = DEG2RAD*get_ddb_angle(vsitetop, resname, "CE1", "ND1", "HD1");
+ a_CE1_ND1_HD1 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CE1", "ND1", "HD1");
}
if (ats[atHE2] != NOTSET)
{
b_NE2_HE2 = get_ddb_bond(vsitetop, resname, "NE2", "HE2");
- a_CE1_NE2_HE2 = DEG2RAD*get_ddb_angle(vsitetop, resname, "CE1", "NE2", "HE2");
+ a_CE1_NE2_HE2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "CE1", "NE2", "HE2");
}
if (ats[atHD2] != NOTSET)
{
b_CD2_HD2 = get_ddb_bond(vsitetop, resname, "CD2", "HD2");
- a_NE2_CD2_HD2 = DEG2RAD*get_ddb_angle(vsitetop, resname, "NE2", "CD2", "HD2");
+ a_NE2_CD2_HD2 = DEG2RAD * get_ddb_angle(vsitetop, resname, "NE2", "CD2", "HD2");
}
if (ats[atHE1] != NOTSET)
{
b_CE1_HE1 = get_ddb_bond(vsitetop, resname, "CE1", "HE1");
- a_NE2_CE1_HE1 = DEG2RAD*get_ddb_angle(vsitetop, resname, "NE2", "CE1", "HE1");
+ a_NE2_CE1_HE1 = DEG2RAD * get_ddb_angle(vsitetop, resname, "NE2", "CE1", "HE1");
}
/* constraints between CG, CE1 and NE1 */
- dCGCE1 = std::sqrt( cosrule(b_CG_ND1, b_ND1_CE1, a_CG_ND1_CE1) );
- dCGNE2 = std::sqrt( cosrule(b_CG_CD2, b_CD2_NE2, a_CG_CD2_NE2) );
+ dCGCE1 = std::sqrt(cosrule(b_CG_ND1, b_ND1_CE1, a_CG_ND1_CE1));
+ dCGNE2 = std::sqrt(cosrule(b_CG_CD2, b_CD2_NE2, a_CG_CD2_NE2));
my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE1);
my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atNE2], dCGNE2);
*/
cosalpha = acosrule(dCGNE2, dCGCE1, b_CE1_NE2);
x[atCE1] = 0;
- y[atCE1] = cosalpha*dCGCE1;
+ y[atCE1] = cosalpha * dCGCE1;
x[atNE2] = 0;
- y[atNE2] = y[atCE1]-b_CE1_NE2;
- sinalpha = std::sqrt(1-cosalpha*cosalpha);
- x[atCG] = -sinalpha*dCGCE1;
+ y[atNE2] = y[atCE1] - b_CE1_NE2;
+ sinalpha = std::sqrt(1 - cosalpha * cosalpha);
+ x[atCG] = -sinalpha * dCGCE1;
y[atCG] = 0;
x[atHE1] = x[atHE2] = x[atHD1] = x[atHD2] = 0;
y[atHE1] = y[atHE2] = y[atHD1] = y[atHD2] = 0;
/* calculate ND1 and CD2 positions from CE1 and NE2 */
- x[atND1] = -b_ND1_CE1*std::sin(a_ND1_CE1_NE2);
- y[atND1] = y[atCE1]-b_ND1_CE1*std::cos(a_ND1_CE1_NE2);
+ x[atND1] = -b_ND1_CE1 * std::sin(a_ND1_CE1_NE2);
+ y[atND1] = y[atCE1] - b_ND1_CE1 * std::cos(a_ND1_CE1_NE2);
- x[atCD2] = -b_CD2_NE2*std::sin(a_CE1_NE2_CD2);
- y[atCD2] = y[atNE2]+b_CD2_NE2*std::cos(a_CE1_NE2_CD2);
+ x[atCD2] = -b_CD2_NE2 * std::sin(a_CE1_NE2_CD2);
+ y[atCD2] = y[atNE2] + b_CD2_NE2 * std::cos(a_CE1_NE2_CD2);
/* And finally the hydrogen positions */
if (ats[atHE1] != NOTSET)
{
- x[atHE1] = x[atCE1] + b_CE1_HE1*std::sin(a_NE2_CE1_HE1);
- y[atHE1] = y[atCE1] - b_CE1_HE1*std::cos(a_NE2_CE1_HE1);
+ x[atHE1] = x[atCE1] + b_CE1_HE1 * std::sin(a_NE2_CE1_HE1);
+ y[atHE1] = y[atCE1] - b_CE1_HE1 * std::cos(a_NE2_CE1_HE1);
}
/* HD2 - first get (ccw) angle from (positive) y-axis */
if (ats[atHD2] != NOTSET)
{
alpha = a_CE1_NE2_CD2 + M_PI - a_NE2_CD2_HD2;
- x[atHD2] = x[atCD2] - b_CD2_HD2*std::sin(alpha);
- y[atHD2] = y[atCD2] + b_CD2_HD2*std::cos(alpha);
+ x[atHD2] = x[atCD2] - b_CD2_HD2 * std::sin(alpha);
+ y[atHD2] = y[atCD2] + b_CD2_HD2 * std::cos(alpha);
}
if (ats[atHD1] != NOTSET)
{
/* HD1 - first get (cw) angle from (positive) y-axis */
alpha = a_ND1_CE1_NE2 + M_PI - a_CE1_ND1_HD1;
- x[atHD1] = x[atND1] - b_ND1_HD1*std::sin(alpha);
- y[atHD1] = y[atND1] - b_ND1_HD1*std::cos(alpha);
+ x[atHD1] = x[atND1] - b_ND1_HD1 * std::sin(alpha);
+ y[atHD1] = y[atND1] - b_ND1_HD1 * std::cos(alpha);
}
if (ats[atHE2] != NOTSET)
{
- x[atHE2] = x[atNE2] + b_NE2_HE2*std::sin(a_CE1_NE2_HE2);
- y[atHE2] = y[atNE2] + b_NE2_HE2*std::cos(a_CE1_NE2_HE2);
+ x[atHE2] = x[atNE2] + b_NE2_HE2 * std::sin(a_CE1_NE2_HE2);
+ y[atHE2] = y[atNE2] + b_NE2_HE2 * std::cos(a_CE1_NE2_HE2);
}
/* Have all coordinates now */
/* calc center-of-mass; keep atoms CG, CE1, NE2 and
* set the rest to vsite3
*/
- mtot = xcom = ycom = 0;
- nvsite = 0;
+ mtot = xcom = ycom = 0;
+ nvsite = 0;
for (i = 0; i < atNR; i++)
{
if (ats[i] != NOTSET)
{
mtot += at->atom[ats[i]].m;
- xcom += x[i]*at->atom[ats[i]].m;
- ycom += y[i]*at->atom[ats[i]].m;
+ xcom += x[i] * at->atom[ats[i]].m;
+ ycom += y[i] * at->atom[ats[i]].m;
if (i != atCG && i != atCE1 && i != atNE2)
{
- at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
- (*vsite_type)[ats[i]] = F_VSITE3;
+ at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
+ (*vsite_type)[ats[i]] = F_VSITE3;
nvsite++;
}
}
}
- if (nvsite+3 != nrfound)
+ if (nvsite + 3 != nrfound)
{
gmx_incons("Generating vsites for HIS");
}
ycom /= mtot;
/* distribute mass so that com stays the same */
- mG = xcom*mtot/x[atCG];
- mrest = mtot-mG;
- mCE1 = (ycom-y[atNE2])*mrest/(y[atCE1]-y[atNE2]);
- mNE2 = mrest-mCE1;
+ mG = xcom * mtot / x[atCG];
+ mrest = mtot - mG;
+ mCE1 = (ycom - y[atNE2]) * mrest / (y[atCE1] - y[atNE2]);
+ mNE2 = mrest - mCE1;
- at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = mG;
+ at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = mG;
at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB = mCE1;
at->atom[ats[atNE2]].m = at->atom[ats[atNE2]].mB = mNE2;
/* HE1 */
if (ats[atHE1] != NOTSET)
{
- calc_vsite3_param(x[atHE1], y[atHE1], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
- x[atCG], y[atCG], &a, &b);
- add_vsite3_param(&plist[F_VSITE3],
- ats[atHE1], ats[atCE1], ats[atNE2], ats[atCG], a, b);
+ calc_vsite3_param(x[atHE1], y[atHE1], x[atCE1], y[atCE1], x[atNE2], y[atNE2], x[atCG],
+ y[atCG], &a, &b);
+ add_vsite3_param(&plist[F_VSITE3], ats[atHE1], ats[atCE1], ats[atNE2], ats[atCG], a, b);
}
/* HE2 */
if (ats[atHE2] != NOTSET)
{
- calc_vsite3_param(x[atHE2], y[atHE2], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
- x[atCG], y[atCG], &a, &b);
- add_vsite3_param(&plist[F_VSITE3],
- ats[atHE2], ats[atNE2], ats[atCE1], ats[atCG], a, b);
+ calc_vsite3_param(x[atHE2], y[atHE2], x[atNE2], y[atNE2], x[atCE1], y[atCE1], x[atCG],
+ y[atCG], &a, &b);
+ add_vsite3_param(&plist[F_VSITE3], ats[atHE2], ats[atNE2], ats[atCE1], ats[atCG], a, b);
}
/* ND1 */
- calc_vsite3_param(x[atND1], y[atND1], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
- x[atCG], y[atCG], &a, &b);
- add_vsite3_param(&plist[F_VSITE3],
- ats[atND1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
+ calc_vsite3_param(x[atND1], y[atND1], x[atNE2], y[atNE2], x[atCE1], y[atCE1], x[atCG], y[atCG],
+ &a, &b);
+ add_vsite3_param(&plist[F_VSITE3], ats[atND1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
/* CD2 */
- calc_vsite3_param(x[atCD2], y[atCD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
- x[atCG], y[atCG], &a, &b);
- add_vsite3_param(&plist[F_VSITE3],
- ats[atCD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
+ calc_vsite3_param(x[atCD2], y[atCD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2], x[atCG], y[atCG],
+ &a, &b);
+ add_vsite3_param(&plist[F_VSITE3], ats[atCD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
/* HD1 */
if (ats[atHD1] != NOTSET)
{
- calc_vsite3_param(x[atHD1], y[atHD1], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
- x[atCG], y[atCG], &a, &b);
- add_vsite3_param(&plist[F_VSITE3],
- ats[atHD1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
+ calc_vsite3_param(x[atHD1], y[atHD1], x[atNE2], y[atNE2], x[atCE1], y[atCE1], x[atCG],
+ y[atCG], &a, &b);
+ add_vsite3_param(&plist[F_VSITE3], ats[atHD1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
}
/* HD2 */
if (ats[atHD2] != NOTSET)
{
- calc_vsite3_param(x[atHD2], y[atHD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
- x[atCG], y[atCG], &a, &b);
- add_vsite3_param(&plist[F_VSITE3],
- ats[atHD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
+ calc_vsite3_param(x[atHD2], y[atHD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2], x[atCG],
+ y[atCG], &a, &b);
+ add_vsite3_param(&plist[F_VSITE3], ats[atHD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
}
return nvsite;
}
{
return FALSE;
}
- switch (abs(vsite_type) )
+ switch (abs(vsite_type))
{
case F_VSITE3:
case F_VSITE3FD:
case F_VSITE3OUT:
case F_VSITE3FAD:
case F_VSITE4FD:
- case F_VSITE4FDN:
- return TRUE;
- default:
- return FALSE;
+ case F_VSITE4FDN: return TRUE;
+ default: return FALSE;
}
}
static char atomnamesuffix[] = "1234";
-void do_vsites(gmx::ArrayRef<const PreprocessResidue> rtpFFDB, PreprocessingAtomTypes *atype,
- t_atoms *at, t_symtab *symtab,
- std::vector<gmx::RVec> *x,
- gmx::ArrayRef<InteractionsOfType> plist, int *vsite_type[], int *cgnr[],
- real mHmult, bool bVsiteAromatics,
- const char *ffdir)
+void do_vsites(gmx::ArrayRef<const PreprocessResidue> rtpFFDB,
+ PreprocessingAtomTypes* atype,
+ t_atoms* at,
+ t_symtab* symtab,
+ std::vector<gmx::RVec>* x,
+ gmx::ArrayRef<InteractionsOfType> plist,
+ int* vsite_type[],
+ int* cgnr[],
+ real mHmult,
+ bool bVsiteAromatics,
+ const char* ffdir)
{
#define MAXATOMSPERRESIDUE 16
- int k, m, i0, ni0, whatres, add_shift, nvsite, nadd;
- int ai, aj, ak, al;
- int nrfound = 0, needed, nrbonds, nrHatoms, Heavy, nrheavies, tpM, tpHeavy;
- int Hatoms[4], heavies[4];
- bool bWARNING, bAddVsiteParam, bFirstWater;
- matrix tmpmat;
- real mHtot, mtot, fact, fact2;
- rvec rpar, rperp, temp;
- char tpname[32], nexttpname[32];
- int *o2n, *newvsite_type, *newcgnr, ats[MAXATOMSPERRESIDUE];
- t_atom *newatom;
- char ***newatomname;
- int cmplength;
- bool isN, planarN, bFound;
+ int k, m, i0, ni0, whatres, add_shift, nvsite, nadd;
+ int ai, aj, ak, al;
+ int nrfound = 0, needed, nrbonds, nrHatoms, Heavy, nrheavies, tpM, tpHeavy;
+ int Hatoms[4], heavies[4];
+ bool bWARNING, bAddVsiteParam, bFirstWater;
+ matrix tmpmat;
+ real mHtot, mtot, fact, fact2;
+ rvec rpar, rperp, temp;
+ char tpname[32], nexttpname[32];
+ int * o2n, *newvsite_type, *newcgnr, ats[MAXATOMSPERRESIDUE];
+ t_atom* newatom;
+ char*** newatomname;
+ int cmplength;
+ bool isN, planarN, bFound;
/* if bVsiteAromatics=TRUE do_vsites will specifically convert atoms in
PHE, TRP, TYR and HIS to a construction of virtual sites */
- enum {
- resPHE, resTRP, resTYR, resHIS, resNR
+ enum
+ {
+ resPHE,
+ resTRP,
+ resTYR,
+ resHIS,
+ resNR
};
- const char *resnms[resNR] = { "PHE", "TRP", "TYR", "HIS" };
+ const char* resnms[resNR] = { "PHE", "TRP", "TYR", "HIS" };
/* Amber03 alternative names for termini */
- const char *resnmsN[resNR] = { "NPHE", "NTRP", "NTYR", "NHIS" };
- const char *resnmsC[resNR] = { "CPHE", "CTRP", "CTYR", "CHIS" };
+ const char* resnmsN[resNR] = { "NPHE", "NTRP", "NTYR", "NHIS" };
+ const char* resnmsC[resNR] = { "CPHE", "CTRP", "CTYR", "CHIS" };
/* HIS can be known as HISH, HIS1, HISA, HID, HIE, HIP, etc. too */
- bool bPartial[resNR] = { FALSE, FALSE, FALSE, TRUE };
+ bool bPartial[resNR] = { FALSE, FALSE, FALSE, TRUE };
/* the atnms for every residue MUST correspond to the enums in the
gen_vsites_* (one for each residue) routines! */
/* also the atom names in atnms MUST be in the same order as in the .rtp! */
- const char *atnms[resNR][MAXATOMSPERRESIDUE+1] = {
+ const char* atnms[resNR][MAXATOMSPERRESIDUE + 1] = {
{ "CG", /* PHE */
- "CD1", "HD1", "CD2", "HD2",
- "CE1", "HE1", "CE2", "HE2",
- "CZ", "HZ", nullptr },
+ "CD1", "HD1", "CD2", "HD2", "CE1", "HE1", "CE2", "HE2", "CZ", "HZ", nullptr },
{ "CB", /* TRP */
- "CG",
- "CD1", "HD1", "CD2",
- "NE1", "HE1", "CE2", "CE3", "HE3",
- "CZ2", "HZ2", "CZ3", "HZ3",
+ "CG", "CD1", "HD1", "CD2", "NE1", "HE1", "CE2", "CE3", "HE3", "CZ2", "HZ2", "CZ3", "HZ3",
"CH2", "HH2", nullptr },
{ "CG", /* TYR */
- "CD1", "HD1", "CD2", "HD2",
- "CE1", "HE1", "CE2", "HE2",
- "CZ", "OH", "HH", nullptr },
+ "CD1", "HD1", "CD2", "HD2", "CE1", "HE1", "CE2", "HE2", "CZ", "OH", "HH", nullptr },
{ "CG", /* HIS */
- "ND1", "HD1", "CD2", "HD2",
- "CE1", "HE1", "NE2", "HE2", nullptr }
+ "ND1", "HD1", "CD2", "HD2", "CE1", "HE1", "NE2", "HE2", nullptr }
};
if (debug)
* force field (rings can be strained).
*/
std::vector<VirtualSiteTopology> vsitetop;
- for (const auto &filename : db)
+ for (const auto& filename : db)
{
read_vsite_database(filename.c_str(), &vsiteconflist, &vsitetop);
}
nvsite = 0;
nadd = 0;
/* we need a marker for which atoms should *not* be renumbered afterwards */
- add_shift = 10*at->nr;
+ add_shift = 10 * at->nr;
/* make arrays where masses can be inserted into */
std::vector<gmx::RVec> newx(at->nr);
snew(newatom, at->nr);
/* make index to tell which residues were already processed */
std::vector<bool> bResProcessed(at->nres);
- ResidueType rt;
+ ResidueType rt;
/* generate vsite constructions */
/* loop over all atoms */
{
resind = at->atom[i].resind;
}
- const char *resnm = *(at->resinfo[resind].name);
+ const char* resnm = *(at->resinfo[resind].name);
/* first check for aromatics to virtualize */
/* don't waste our effort on DNA, water etc. */
/* Only do the vsite aromatic stuff when we reach the
* CA atom, since there might be an X2/X3 group on the
* N-terminus that must be treated first.
*/
- if (bVsiteAromatics &&
- (strcmp(*(at->atomname[i]), "CA") == 0) &&
- !bResProcessed[resind] &&
- rt.namedResidueHasType(*(at->resinfo[resind].name), "Protein") )
+ if (bVsiteAromatics && (strcmp(*(at->atomname[i]), "CA") == 0) && !bResProcessed[resind]
+ && rt.namedResidueHasType(*(at->resinfo[resind].name), "Protein"))
{
/* mark this residue */
bResProcessed[resind] = TRUE;
for (int j = 0; j < resNR && whatres == NOTSET; j++)
{
- cmplength = bPartial[j] ? strlen(resnm)-1 : strlen(resnm);
+ cmplength = bPartial[j] ? strlen(resnm) - 1 : strlen(resnm);
- bFound = ((gmx::equalCaseInsensitive(resnm, resnms[j], cmplength)) ||
- (gmx::equalCaseInsensitive(resnm, resnmsN[j], cmplength)) ||
- (gmx::equalCaseInsensitive(resnm, resnmsC[j], cmplength)));
+ bFound = ((gmx::equalCaseInsensitive(resnm, resnms[j], cmplength))
+ || (gmx::equalCaseInsensitive(resnm, resnmsN[j], cmplength))
+ || (gmx::equalCaseInsensitive(resnm, resnmsC[j], cmplength)));
if (bFound)
{
/* now k is number of atom names in atnms[j] */
if (j == resHIS)
{
- needed = k-3;
+ needed = k - 3;
}
else
{
}
if (nrfound < needed)
{
- gmx_fatal(FARGS, "not enough atoms found (%d, need %d) in "
+ gmx_fatal(FARGS,
+ "not enough atoms found (%d, need %d) in "
"residue %s %d while\n "
"generating aromatics virtual site construction",
nrfound, needed, resnm, at->resinfo[resind].nr);
case resPHE:
if (debug)
{
- fprintf(stderr, "PHE at %d\n", o2n[ats[0]]+1);
+ fprintf(stderr, "PHE at %d\n", o2n[ats[0]] + 1);
}
nvsite += gen_vsites_phe(at, vsite_type, plist, nrfound, ats, vsitetop);
break;
case resTRP:
if (debug)
{
- fprintf(stderr, "TRP at %d\n", o2n[ats[0]]+1);
+ fprintf(stderr, "TRP at %d\n", o2n[ats[0]] + 1);
}
nvsite += gen_vsites_trp(atype, &newx, &newatom, &newatomname, &o2n,
- &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
- at, vsite_type, plist, nrfound, ats, add_shift, vsitetop);
+ &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr, at,
+ vsite_type, plist, nrfound, ats, add_shift, vsitetop);
break;
case resTYR:
if (debug)
{
- fprintf(stderr, "TYR at %d\n", o2n[ats[0]]+1);
+ fprintf(stderr, "TYR at %d\n", o2n[ats[0]] + 1);
}
nvsite += gen_vsites_tyr(atype, &newx, &newatom, &newatomname, &o2n,
- &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
- at, vsite_type, plist, nrfound, ats, add_shift, vsitetop);
+ &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr, at,
+ vsite_type, plist, nrfound, ats, add_shift, vsitetop);
break;
case resHIS:
if (debug)
{
- fprintf(stderr, "HIS at %d\n", o2n[ats[0]]+1);
+ fprintf(stderr, "HIS at %d\n", o2n[ats[0]] + 1);
}
nvsite += gen_vsites_his(at, vsite_type, plist, nrfound, ats, vsitetop);
break;
case NOTSET:
/* this means this residue won't be processed */
break;
- default:
- gmx_fatal(FARGS, "DEATH HORROR in do_vsites (%s:%d)",
- __FILE__, __LINE__);
+ default: gmx_fatal(FARGS, "DEATH HORROR in do_vsites (%s:%d)", __FILE__, __LINE__);
} /* switch whatres */
/* skip back to beginning of residue */
- while (i > 0 && at->atom[i-1].resind == resind)
+ while (i > 0 && at->atom[i - 1].resind == resind)
{
i--;
}
/* now process the rest of the hydrogens */
/* only process hydrogen atoms which are not already set */
- if ( ((*vsite_type)[i] == NOTSET) && is_hydrogen(*(at->atomname[i])))
+ if (((*vsite_type)[i] == NOTSET) && is_hydrogen(*(at->atomname[i])))
{
/* find heavy atom, count #bonds from it and #H atoms bound to it
and return H atom numbers (Hatoms) and heavy atom numbers (heavies) */
- count_bonds(i, &plist[F_BONDS], at->atomname,
- &nrbonds, &nrHatoms, Hatoms, &Heavy, &nrheavies, heavies);
+ count_bonds(i, &plist[F_BONDS], at->atomname, &nrbonds, &nrHatoms, Hatoms, &Heavy,
+ &nrheavies, heavies);
/* get Heavy atom type */
tpHeavy = get_atype(Heavy, at, rtpFFDB, &rt);
strcpy(tpname, atype->atomNameFromAtomType(tpHeavy));
{
switch (nrbonds)
{
- case 2: /* -O-H */
- (*vsite_type)[i] = F_BONDS;
- break;
- case 3: /* =CH-, -NH- or =NH+- */
- (*vsite_type)[i] = F_VSITE3FD;
- break;
+ case 2: /* -O-H */ (*vsite_type)[i] = F_BONDS; break;
+ case 3: /* =CH-, -NH- or =NH+- */ (*vsite_type)[i] = F_VSITE3FD; break;
case 4: /* --CH- (tert) */
/* The old type 4FD had stability issues, so
* all new constructs should use 4FDN
}
break;
- default: /* nrbonds != 2, 3 or 4 */
- bWARNING = TRUE;
+ default: /* nrbonds != 2, 3 or 4 */ bWARNING = TRUE;
}
-
}
- else if ( (nrHatoms == 2) && (nrbonds == 2) &&
- (at->atom[Heavy].atomnumber == 8) )
+ else if ((nrHatoms == 2) && (nrbonds == 2) && (at->atom[Heavy].atomnumber == 8))
{
bAddVsiteParam = FALSE; /* this is water: skip these hydrogens */
if (bFirstWater)
bFirstWater = FALSE;
if (debug)
{
- fprintf(debug,
- "Not converting hydrogens in water to virtual sites\n");
+ fprintf(debug, "Not converting hydrogens in water to virtual sites\n");
}
}
}
- else if ( (nrHatoms == 2) && (nrbonds == 4) )
+ else if ((nrHatoms == 2) && (nrbonds == 4))
{
/* -CH2- , -NH2+- */
(*vsite_type)[Hatoms[0]] = F_VSITE3OUT;
isN = planarN = FALSE;
if ((nrHatoms == 2) && ((*at->atomname[Heavy])[0] == 'N'))
{
- isN = TRUE;
- int j = nitrogen_is_planar(vsiteconflist, tpname);
+ isN = TRUE;
+ int j = nitrogen_is_planar(vsiteconflist, tpname);
if (j < 0)
{
gmx_fatal(FARGS, "No vsite database NH2 entry for type %s\n", tpname);
}
planarN = (j == 1);
}
- if ( (nrHatoms == 2) && (nrbonds == 3) && ( !isN || planarN ) )
+ if ((nrHatoms == 2) && (nrbonds == 3) && (!isN || planarN))
{
/* =CH2 or, if it is a nitrogen NH2, it is a planar one */
(*vsite_type)[Hatoms[0]] = F_VSITE3FAD;
(*vsite_type)[Hatoms[1]] = -F_VSITE3FAD;
}
- else if ( ( (nrHatoms == 2) && (nrbonds == 3) &&
- ( isN && !planarN ) ) ||
- ( (nrHatoms == 3) && (nrbonds == 4) ) )
+ else if (((nrHatoms == 2) && (nrbonds == 3) && (isN && !planarN))
+ || ((nrHatoms == 3) && (nrbonds == 4)))
{
/* CH3, NH3 or non-planar NH2 group */
- int Hat_vsite_type[3] = { F_VSITE3, F_VSITE3OUT, F_VSITE3OUT };
- bool Hat_SwapParity[3] = { FALSE, TRUE, FALSE };
+ int Hat_vsite_type[3] = { F_VSITE3, F_VSITE3OUT, F_VSITE3OUT };
+ bool Hat_SwapParity[3] = { FALSE, TRUE, FALSE };
if (debug)
{
- fprintf(stderr, "-XH3 or nonplanar NH2 group at %d\n", i+1);
+ fprintf(stderr, "-XH3 or nonplanar NH2 group at %d\n", i + 1);
}
bAddVsiteParam = FALSE; /* we'll do this ourselves! */
/* -NH2 (umbrella), -NH3+ or -CH3 */
- (*vsite_type)[Heavy] = F_VSITE3;
+ (*vsite_type)[Heavy] = F_VSITE3;
for (int j = 0; j < nrHatoms; j++)
{
(*vsite_type)[Hatoms[j]] = Hat_vsite_type[j];
}
/* get dummy mass type from first char of heavy atom type (N or C) */
- strcpy(nexttpname, atype->atomNameFromAtomType(get_atype(heavies[0], at, rtpFFDB, &rt)));
+ strcpy(nexttpname,
+ atype->atomNameFromAtomType(get_atype(heavies[0], at, rtpFFDB, &rt)));
std::string ch = get_dummymass_name(vsiteconflist, tpname, nexttpname);
std::string name;
if (ch.empty())
{
if (!db.empty())
{
- 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);
+ 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);
}
else
{
- 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);
+ 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);
}
}
else
/* make space for 2 masses: shift all atoms starting with 'Heavy' */
#define NMASS 2
i0 = Heavy;
- ni0 = i0+nadd;
+ ni0 = i0 + nadd;
if (debug)
{
- fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, o2n[i0]+1);
+ fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, o2n[i0] + 1);
}
nadd += NMASS;
for (int j = i0; j < at->nr; j++)
{
- o2n[j] = j+nadd;
+ o2n[j] = j + nadd;
}
- newx.resize(at->nr+nadd);
- srenew(newatom, at->nr+nadd);
- srenew(newatomname, at->nr+nadd);
- srenew(newvsite_type, at->nr+nadd);
- srenew(newcgnr, at->nr+nadd);
+ newx.resize(at->nr + nadd);
+ srenew(newatom, at->nr + nadd);
+ srenew(newatomname, at->nr + nadd);
+ srenew(newvsite_type, at->nr + nadd);
+ srenew(newcgnr, at->nr + nadd);
for (int j = 0; j < NMASS; j++)
{
- newatomname[at->nr+nadd-1-j] = nullptr;
+ newatomname[at->nr + nadd - 1 - j] = nullptr;
}
/* calculate starting position for the masses */
/* get atom masses, and set Heavy and Hatoms mass to zero */
for (int j = 0; j < nrHatoms; j++)
{
- mHtot += get_amass(Hatoms[j], at, rtpFFDB, &rt);
+ mHtot += get_amass(Hatoms[j], at, rtpFFDB, &rt);
at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
}
mtot = mHtot + get_amass(Heavy, at, rtpFFDB, &rt);
{
mHtot *= mHmult;
}
- fact2 = mHtot/mtot;
+ fact2 = mHtot / mtot;
fact = std::sqrt(fact2);
/* generate vectors parallel and perpendicular to rotational axis:
* rpar = Heavy -> Hcom
{
rvec_inc(rpar, (*x)[Hatoms[j]]);
}
- svmul(1.0/nrHatoms, rpar, rpar); /* rpar = ( H1+H2+H3 ) / 3 */
- rvec_dec(rpar, (*x)[Heavy]); /* - Heavy */
+ svmul(1.0 / nrHatoms, rpar, rpar); /* rpar = ( H1+H2+H3 ) / 3 */
+ rvec_dec(rpar, (*x)[Heavy]); /* - Heavy */
rvec_sub((*x)[Hatoms[0]], (*x)[Heavy], rperp);
- rvec_dec(rperp, rpar); /* rperp = H1 - Heavy - rpar */
+ rvec_dec(rperp, rpar); /* rperp = H1 - Heavy - rpar */
/* calc mass positions */
svmul(fact2, rpar, temp);
for (int j = 0; (j < NMASS); j++) /* xM = xN + fact2 * rpar +/- fact * rperp */
{
- rvec_add((*x)[Heavy], temp, newx[ni0+j]);
+ rvec_add((*x)[Heavy], temp, newx[ni0 + j]);
}
svmul(fact, rperp, temp);
- rvec_inc(newx[ni0 ], temp);
- rvec_dec(newx[ni0+1], temp);
+ rvec_inc(newx[ni0], temp);
+ rvec_dec(newx[ni0 + 1], temp);
/* set atom parameters for the masses */
for (int j = 0; (j < NMASS); j++)
{
/* make name: "M??#" or "M?#" (? is atomname, # is number) */
name[0] = 'M';
int k;
- for (k = 0; (*at->atomname[Heavy])[k] && ( k < NMASS ); k++)
+ for (k = 0; (*at->atomname[Heavy])[k] && (k < NMASS); k++)
{
- name[k+1] = (*at->atomname[Heavy])[k];
+ name[k + 1] = (*at->atomname[Heavy])[k];
}
- name[k+1] = atomnamesuffix[j];
- name[k+2] = '\0';
- newatomname[ni0+j] = put_symtab(symtab, name.c_str());
- newatom[ni0+j].m = newatom[ni0+j].mB = mtot/NMASS;
- newatom[ni0+j].q = newatom[ni0+j].qB = 0.0;
- newatom[ni0+j].type = newatom[ni0+j].typeB = tpM;
- newatom[ni0+j].ptype = eptAtom;
- newatom[ni0+j].resind = at->atom[i0].resind;
- newatom[ni0+j].elem[0] = 'M';
- newatom[ni0+j].elem[1] = '\0';
- newvsite_type[ni0+j] = NOTSET;
- newcgnr[ni0+j] = (*cgnr)[i0];
+ name[k + 1] = atomnamesuffix[j];
+ name[k + 2] = '\0';
+ newatomname[ni0 + j] = put_symtab(symtab, name.c_str());
+ newatom[ni0 + j].m = newatom[ni0 + j].mB = mtot / NMASS;
+ newatom[ni0 + j].q = newatom[ni0 + j].qB = 0.0;
+ newatom[ni0 + j].type = newatom[ni0 + j].typeB = tpM;
+ newatom[ni0 + j].ptype = eptAtom;
+ newatom[ni0 + j].resind = at->atom[i0].resind;
+ newatom[ni0 + j].elem[0] = 'M';
+ newatom[ni0 + j].elem[1] = '\0';
+ newvsite_type[ni0 + j] = NOTSET;
+ newcgnr[ni0 + j] = (*cgnr)[i0];
}
/* add constraints between dummy masses and to heavies[0] */
/* 'add_shift' says which atoms won't be renumbered afterwards */
- my_add_param(&(plist[F_CONSTRNC]), heavies[0], add_shift+ni0, NOTSET);
- my_add_param(&(plist[F_CONSTRNC]), heavies[0], add_shift+ni0+1, NOTSET);
- my_add_param(&(plist[F_CONSTRNC]), add_shift+ni0, add_shift+ni0+1, NOTSET);
+ my_add_param(&(plist[F_CONSTRNC]), heavies[0], add_shift + ni0, NOTSET);
+ my_add_param(&(plist[F_CONSTRNC]), heavies[0], add_shift + ni0 + 1, NOTSET);
+ my_add_param(&(plist[F_CONSTRNC]), add_shift + ni0, add_shift + ni0 + 1, NOTSET);
/* generate Heavy, H1, H2 and H3 from M1, M2 and heavies[0] */
/* note that vsite_type cannot be NOTSET, because we just set it */
- add_vsite3_atoms (&plist[(*vsite_type)[Heavy]],
- Heavy, heavies[0], add_shift+ni0, add_shift+ni0+1,
- FALSE);
+ add_vsite3_atoms(&plist[(*vsite_type)[Heavy]], Heavy, heavies[0],
+ add_shift + ni0, add_shift + ni0 + 1, FALSE);
for (int j = 0; j < nrHatoms; j++)
{
- add_vsite3_atoms(&plist[(*vsite_type)[Hatoms[j]]],
- Hatoms[j], heavies[0], add_shift+ni0, add_shift+ni0+1,
- Hat_SwapParity[j]);
+ add_vsite3_atoms(&plist[(*vsite_type)[Hatoms[j]]], Hatoms[j], heavies[0],
+ add_shift + ni0, add_shift + ni0 + 1, Hat_SwapParity[j]);
}
#undef NMASS
}
{
bWARNING = TRUE;
}
-
}
if (bWARNING)
{
- gmx_fatal(FARGS, "Cannot convert atom %d %s (bound to a heavy atom "
+ gmx_fatal(FARGS,
+ "Cannot convert atom %d %s (bound to a heavy atom "
"%s with \n"
" %d bonds and %d bound hydrogens atoms) to virtual site\n",
- i+1, *(at->atomname[i]), tpname, nrbonds, nrHatoms);
+ i + 1, *(at->atomname[i]), tpname, nrbonds, nrHatoms);
}
if (bAddVsiteParam)
{
/* add vsite parameters to topology,
also get rid of negative vsite_types */
- add_vsites(plist, (*vsite_type), Heavy, nrHatoms, Hatoms,
- nrheavies, heavies);
+ add_vsites(plist, (*vsite_type), Heavy, nrHatoms, Hatoms, nrheavies, heavies);
/* transfer mass of virtual site to Heavy atom */
for (int j = 0; j < nrHatoms; j++)
{
if (is_vsite((*vsite_type)[Hatoms[j]]))
{
- at->atom[Heavy].m += at->atom[Hatoms[j]].m;
+ at->atom[Heavy].m += at->atom[Hatoms[j]].m;
at->atom[Heavy].mB = at->atom[Heavy].m;
at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
}
nvsite += nrHatoms;
if (debug)
{
- fprintf(debug, "atom %d: ", o2n[i]+1);
+ fprintf(debug, "atom %d: ", o2n[i] + 1);
print_bonds(debug, o2n, nrHatoms, Hatoms, Heavy, nrheavies, heavies);
}
} /* if vsite NOTSET & is hydrogen */
- } /* for i < at->nr */
+ } /* for i < at->nr */
if (debug)
{
fprintf(debug, "Before inserting new atoms:\n");
for (int i = 0; i < at->nr; i++)
{
- fprintf(debug, "%4d %4d %4s %4d %4s %6d %-10s\n", i+1, o2n[i]+1,
- at->atomname[i] ? *(at->atomname[i]) : "(NULL)",
- at->resinfo[at->atom[i].resind].nr,
- at->resinfo[at->atom[i].resind].name ?
- *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
+ fprintf(debug, "%4d %4d %4s %4d %4s %6d %-10s\n", i + 1, o2n[i] + 1,
+ at->atomname[i] ? *(at->atomname[i]) : "(NULL)", at->resinfo[at->atom[i].resind].nr,
+ at->resinfo[at->atom[i].resind].name ? *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
(*cgnr)[i],
- ((*vsite_type)[i] == NOTSET) ?
- "NOTSET" : interaction_function[(*vsite_type)[i]].name);
+ ((*vsite_type)[i] == NOTSET) ? "NOTSET" : interaction_function[(*vsite_type)[i]].name);
}
fprintf(debug, "new atoms to be inserted:\n");
- for (int i = 0; i < at->nr+nadd; i++)
+ for (int i = 0; i < at->nr + nadd; i++)
{
if (newatomname[i])
{
- fprintf(debug, "%4d %4s %4d %6d %-10s\n", i+1,
- newatomname[i] ? *(newatomname[i]) : "(NULL)",
- newatom[i].resind, newcgnr[i],
- (newvsite_type[i] == NOTSET) ?
- "NOTSET" : interaction_function[newvsite_type[i]].name);
+ fprintf(debug, "%4d %4s %4d %6d %-10s\n", i + 1,
+ newatomname[i] ? *(newatomname[i]) : "(NULL)", newatom[i].resind, newcgnr[i],
+ (newvsite_type[i] == NOTSET) ? "NOTSET"
+ : interaction_function[newvsite_type[i]].name);
}
}
}
/* add all original atoms to the new arrays, using o2n index array */
for (int i = 0; i < at->nr; i++)
{
- newatomname [o2n[i]] = at->atomname [i];
- newatom [o2n[i]] = at->atom [i];
+ newatomname[o2n[i]] = at->atomname[i];
+ newatom[o2n[i]] = at->atom[i];
newvsite_type[o2n[i]] = (*vsite_type)[i];
- newcgnr [o2n[i]] = (*cgnr) [i];
+ newcgnr[o2n[i]] = (*cgnr)[i];
copy_rvec((*x)[i], newx[o2n[i]]);
}
/* throw away old atoms */
sfree(*vsite_type);
sfree(*cgnr);
/* put in the new ones */
- at->nr += nadd;
+ at->nr += nadd;
at->atom = newatom;
at->atomname = newatomname;
*vsite_type = newvsite_type;
*x = newx;
if (at->nr > add_shift)
{
- gmx_fatal(FARGS, "Added impossible amount of dummy masses "
- "(%d on a total of %d atoms)\n", nadd, at->nr-nadd);
+ gmx_fatal(FARGS,
+ "Added impossible amount of dummy masses "
+ "(%d on a total of %d atoms)\n",
+ nadd, at->nr - nadd);
}
if (debug)
fprintf(debug, "After inserting new atoms:\n");
for (int i = 0; i < at->nr; i++)
{
- fprintf(debug, "%4d %4s %4d %4s %6d %-10s\n", i+1,
- at->atomname[i] ? *(at->atomname[i]) : "(NULL)",
- at->resinfo[at->atom[i].resind].nr,
- at->resinfo[at->atom[i].resind].name ?
- *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
+ fprintf(debug, "%4d %4s %4d %4s %6d %-10s\n", i + 1,
+ at->atomname[i] ? *(at->atomname[i]) : "(NULL)", at->resinfo[at->atom[i].resind].nr,
+ at->resinfo[at->atom[i].resind].name ? *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
(*cgnr)[i],
- ((*vsite_type)[i] == NOTSET) ?
- "NOTSET" : interaction_function[(*vsite_type)[i]].name);
+ ((*vsite_type)[i] == NOTSET) ? "NOTSET" : interaction_function[(*vsite_type)[i]].name);
}
}
/* now renumber all the interactions because of the added atoms */
for (int ftype = 0; ftype < F_NRE; ftype++)
{
- InteractionsOfType *params = &(plist[ftype]);
+ InteractionsOfType* params = &(plist[ftype]);
if (debug)
{
- fprintf(debug, "Renumbering %zu %s\n", params->size(),
- interaction_function[ftype].longname);
+ fprintf(debug, "Renumbering %zu %s\n", params->size(), interaction_function[ftype].longname);
}
/* Horrible hacks needed here to get this to work */
for (auto parm = params->interactionTypes.begin(); parm != params->interactionTypes.end(); parm++)
{
gmx::ArrayRef<const int> atomNumbers(parm->atoms());
- std::vector<int> newAtomNumber;
+ std::vector<int> newAtomNumber;
for (int j = 0; j < NRAL(ftype); j++)
{
if (atomNumbers[j] >= add_shift)
{
if (debug)
{
- fprintf(debug, " [%d -> %d]", atomNumbers[j],
- atomNumbers[j]-add_shift);
+ fprintf(debug, " [%d -> %d]", atomNumbers[j], atomNumbers[j] - add_shift);
}
- newAtomNumber.emplace_back(atomNumbers[j]-add_shift);
+ newAtomNumber.emplace_back(atomNumbers[j] - add_shift);
}
else
{
if (debug)
{
- fprintf(debug, " [%d -> %d]", atomNumbers[j],
- o2n[atomNumbers[j]]);
+ fprintf(debug, " [%d -> %d]", atomNumbers[j], o2n[atomNumbers[j]]);
}
newAtomNumber.emplace_back(o2n[atomNumbers[j]]);
}
}
}
/* sort constraint parameters */
- InteractionsOfType *params = &(plist[F_CONSTRNC]);
- for (auto &type : params->interactionTypes)
+ InteractionsOfType* params = &(plist[F_CONSTRNC]);
+ for (auto& type : params->interactionTypes)
{
type.sortAtomIds();
}
fprintf(stderr, "Added %zu new constraints\n", plist[F_CONSTRNC].size());
}
-void do_h_mass(InteractionsOfType *psb, int vsite_type[], t_atoms *at, real mHmult,
- bool bDeuterate)
+void do_h_mass(InteractionsOfType* psb, int vsite_type[], t_atoms* at, real mHmult, bool bDeuterate)
{
/* loop over all atoms */
for (int i = 0; i < at->nr; i++)
{
/* adjust masses if i is hydrogen and not a virtual site */
- if (!is_vsite(vsite_type[i]) && is_hydrogen(*(at->atomname[i])) )
+ if (!is_vsite(vsite_type[i]) && is_hydrogen(*(at->atomname[i])))
{
/* find bonded heavy atom */
int a = NOTSET;
- for (auto parm = psb->interactionTypes.begin(); (parm != psb->interactionTypes.end()) && (a == NOTSET); parm++)
+ for (auto parm = psb->interactionTypes.begin();
+ (parm != psb->interactionTypes.end()) && (a == NOTSET); parm++)
{
/* if other atom is not a virtual site, it is the one we want */
- if ( (parm->ai() == i) &&
- !is_vsite(vsite_type[parm->aj()]) )
+ if ((parm->ai() == i) && !is_vsite(vsite_type[parm->aj()]))
{
a = parm->aj();
}
- else if ( (parm->aj() == i) &&
- !is_vsite(vsite_type[parm->ai()]) )
+ else if ((parm->aj() == i) && !is_vsite(vsite_type[parm->ai()]))
{
a = parm->ai();
}
}
if (a == NOTSET)
{
- gmx_fatal(FARGS, "Unbound hydrogen atom (%d) found while adjusting mass",
- i+1);
+ gmx_fatal(FARGS, "Unbound hydrogen atom (%d) found while adjusting mass", i + 1);
}
/* adjust mass of i (hydrogen) with mHmult
and correct mass of a (bonded atom) with same amount */
if (!bDeuterate)
{
- at->atom[a].m -= (mHmult-1.0)*at->atom[i].m;
- at->atom[a].mB -= (mHmult-1.0)*at->atom[i].m;
+ at->atom[a].m -= (mHmult - 1.0) * at->atom[i].m;
+ at->atom[a].mB -= (mHmult - 1.0) * at->atom[i].m;
}
- at->atom[i].m *= mHmult;
+ at->atom[i].m *= mHmult;
at->atom[i].mB *= mHmult;
}
}
biod.pnpi.spb.ru / alexxy / gromacs.git / blobdiff