3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
10 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
11 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
12 * Copyright (c) 2001-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
15 * This program is free software; you can redistribute it and/or
16 * modify it under the terms of the GNU General Public License
17 * as published by the Free Software Foundation; either version 2
18 * of the License, or (at your option) any later version.
20 * If you want to redistribute modifications, please consider that
21 * scientific software is very special. Version control is crucial -
22 * bugs must be traceable. We will be happy to consider code for
23 * inclusion in the official distribution, but derived work must not
24 * be called official GROMACS. Details are found in the README & COPYING
25 * files - if they are missing, get the official version at www.gromacs.org.
27 * To help us fund GROMACS development, we humbly ask that you cite
28 * the papers on the package - you can find them in the top README file.
30 * For more info, check our website at http://www.gromacs.org
33 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
43 #include "gen_vsite.h"
53 #include "gpp_atomtype.h"
54 #include "fflibutil.h"
58 #define OPENDIR '[' /* starting sign for directive */
59 #define CLOSEDIR ']' /* ending sign for directive */
62 char atomtype[MAXNAME]; /* Type for the XH3/XH2 atom */
63 gmx_bool isplanar; /* If true, the atomtype above and the three connected
64 * ones are in a planar geometry. The two next entries
65 * are undefined in that case
67 int nhydrogens; /* number of connected hydrogens */
68 char nextheavytype[MAXNAME]; /* Type for the heavy atom bonded to XH2/XH3 */
69 char dummymass[MAXNAME]; /* The type of MNH* or MCH3* dummy mass to use */
73 /* Structure to represent average bond and angles values in vsite aromatic
74 * residues. Note that these are NOT necessarily the bonds and angles from the
75 * forcefield; many forcefields (like Amber, OPLS) have some inherent strain in
76 * 5-rings (i.e. the sum of angles is !=540, but impropers keep it planar)
79 char resname[MAXNAME];
82 struct vsitetop_bond {
86 } *bond; /* list of bonds */
87 struct vsitetop_angle {
92 } *angle; /* list of angles */
97 DDB_CH3, DDB_NH3, DDB_NH2, DDB_PHE, DDB_TYR,
98 DDB_TRP, DDB_HISA, DDB_HISB, DDB_HISH, DDB_DIR_NR
101 typedef char t_dirname[STRLEN];
103 static const t_dirname ddb_dirnames[DDB_DIR_NR] = {
115 static int ddb_name2dir(char *name)
117 /* Translate a directive name to the number of the directive.
118 * HID/HIE/HIP names are translated to the ones we use in Gromacs.
125 for (i = 0; i < DDB_DIR_NR && index < 0; i++)
127 if (!gmx_strcasecmp(name, ddb_dirnames[i]))
137 static void read_vsite_database(const char *ddbname,
138 t_vsiteconf **pvsiteconflist, int *nvsiteconf,
139 t_vsitetop **pvsitetoplist, int *nvsitetop)
141 /* This routine is a quick hack to fix the problem with hardcoded atomtypes
142 * and aromatic vsite parameters by reading them from a ff???.vsd file.
144 * The file can contain sections [ NH3 ], [ CH3 ], [ NH2 ], and ring residue names.
145 * For the NH3 and CH3 section each line has three fields. The first is the atomtype
146 * (nb: not bonded type) of the N/C atom to be replaced, the second field is
147 * the type of the next heavy atom it is bonded to, and the third field the type
148 * of dummy mass that will be used for this group.
150 * If the NH2 group planar (sp2 N) a different vsite construct is used, so in this
151 * case the second field should just be the word planar.
157 int i, j, n, k, nvsite, ntop, curdir, prevdir;
158 t_vsiteconf *vsiteconflist;
159 t_vsitetop *vsitetoplist;
161 char s1[MAXNAME], s2[MAXNAME], s3[MAXNAME], s4[MAXNAME];
163 ddb = libopen(ddbname);
165 nvsite = *nvsiteconf;
166 vsiteconflist = *pvsiteconflist;
168 vsitetoplist = *pvsitetoplist;
172 snew(vsiteconflist, 1);
173 snew(vsitetoplist, 1);
175 while (fgets2(pline, STRLEN-2, ddb) != NULL)
177 strip_comment(pline);
179 if (strlen(pline) > 0)
181 if (pline[0] == OPENDIR)
183 strncpy(dirstr, pline+1, STRLEN-2);
184 if ((ch = strchr (dirstr, CLOSEDIR)) != NULL)
190 if (!gmx_strcasecmp(dirstr, "HID") ||
191 !gmx_strcasecmp(dirstr, "HISD"))
193 sprintf(dirstr, "HISA");
195 else if (!gmx_strcasecmp(dirstr, "HIE") ||
196 !gmx_strcasecmp(dirstr, "HISE"))
198 sprintf(dirstr, "HISB");
200 else if (!gmx_strcasecmp(dirstr, "HIP"))
202 sprintf(dirstr, "HISH");
205 curdir = ddb_name2dir(dirstr);
208 gmx_fatal(FARGS, "Invalid directive %s in vsite database %s",
217 gmx_fatal(FARGS, "First entry in vsite database must be a directive.\n");
222 n = sscanf(pline, "%s%s%s", s1, s2, s3);
223 if (n < 3 && !gmx_strcasecmp(s2, "planar"))
225 srenew(vsiteconflist, nvsite+1);
226 strncpy(vsiteconflist[nvsite].atomtype, s1, MAXNAME-1);
227 vsiteconflist[nvsite].isplanar = TRUE;
228 vsiteconflist[nvsite].nextheavytype[0] = 0;
229 vsiteconflist[nvsite].dummymass[0] = 0;
230 vsiteconflist[nvsite].nhydrogens = 2;
235 srenew(vsiteconflist, (nvsite+1));
236 strncpy(vsiteconflist[nvsite].atomtype, s1, MAXNAME-1);
237 vsiteconflist[nvsite].isplanar = FALSE;
238 strncpy(vsiteconflist[nvsite].nextheavytype, s2, MAXNAME-1);
239 strncpy(vsiteconflist[nvsite].dummymass, s3, MAXNAME-1);
240 if (curdir == DDB_NH2)
242 vsiteconflist[nvsite].nhydrogens = 2;
246 vsiteconflist[nvsite].nhydrogens = 3;
252 gmx_fatal(FARGS, "Not enough directives in vsite database line: %s\n", pline);
262 while ((i < ntop) && gmx_strcasecmp(dirstr, vsitetoplist[i].resname))
266 /* Allocate a new topology entry if this is a new residue */
269 srenew(vsitetoplist, ntop+1);
270 ntop++; /* i still points to current vsite topology entry */
271 strncpy(vsitetoplist[i].resname, dirstr, MAXNAME-1);
272 vsitetoplist[i].nbonds = vsitetoplist[i].nangles = 0;
273 snew(vsitetoplist[i].bond, 1);
274 snew(vsitetoplist[i].angle, 1);
276 n = sscanf(pline, "%s%s%s%s", s1, s2, s3, s4);
280 k = vsitetoplist[i].nbonds++;
281 srenew(vsitetoplist[i].bond, k+1);
282 strncpy(vsitetoplist[i].bond[k].atom1, s1, MAXNAME-1);
283 strncpy(vsitetoplist[i].bond[k].atom2, s2, MAXNAME-1);
284 vsitetoplist[i].bond[k].value = strtod(s3, NULL);
289 k = vsitetoplist[i].nangles++;
290 srenew(vsitetoplist[i].angle, k+1);
291 strncpy(vsitetoplist[i].angle[k].atom1, s1, MAXNAME-1);
292 strncpy(vsitetoplist[i].angle[k].atom2, s2, MAXNAME-1);
293 strncpy(vsitetoplist[i].angle[k].atom3, s3, MAXNAME-1);
294 vsitetoplist[i].angle[k].value = strtod(s4, NULL);
298 gmx_fatal(FARGS, "Need 3 or 4 values to specify bond/angle values in %s: %s\n", ddbname, pline);
302 gmx_fatal(FARGS, "Didnt find a case for directive %s in read_vsite_database\n", dirstr);
308 *pvsiteconflist = vsiteconflist;
309 *pvsitetoplist = vsitetoplist;
310 *nvsiteconf = nvsite;
316 static int nitrogen_is_planar(t_vsiteconf vsiteconflist[], int nvsiteconf, char atomtype[])
318 /* Return 1 if atomtype exists in database list and is planar, 0 if not,
319 * and -1 if not found.
322 gmx_bool found = FALSE;
323 for (i = 0; i < nvsiteconf && !found; i++)
325 found = (!gmx_strcasecmp(vsiteconflist[i].atomtype, atomtype) && (vsiteconflist[i].nhydrogens == 2));
329 res = (vsiteconflist[i-1].isplanar == TRUE);
339 static char *get_dummymass_name(t_vsiteconf vsiteconflist[], int nvsiteconf, char atom[], char nextheavy[])
341 /* Return the dummy mass name if found, or NULL if not set in ddb database */
343 gmx_bool found = FALSE;
344 for (i = 0; i < nvsiteconf && !found; i++)
346 found = (!gmx_strcasecmp(vsiteconflist[i].atomtype, atom) &&
347 !gmx_strcasecmp(vsiteconflist[i].nextheavytype, nextheavy));
351 return vsiteconflist[i-1].dummymass;
361 static real get_ddb_bond(t_vsitetop *vsitetop, int nvsitetop,
363 const char atom1[], const char atom2[])
368 while (i < nvsitetop && gmx_strcasecmp(res, vsitetop[i].resname))
374 gmx_fatal(FARGS, "No vsite information for residue %s found in vsite database.\n", res);
377 while (j < vsitetop[i].nbonds &&
378 ( strcmp(atom1, vsitetop[i].bond[j].atom1) || strcmp(atom2, vsitetop[i].bond[j].atom2)) &&
379 ( strcmp(atom2, vsitetop[i].bond[j].atom1) || strcmp(atom1, vsitetop[i].bond[j].atom2)))
383 if (j == vsitetop[i].nbonds)
385 gmx_fatal(FARGS, "Couldnt find bond %s-%s for residue %s in vsite database.\n", atom1, atom2, res);
388 return vsitetop[i].bond[j].value;
392 static real get_ddb_angle(t_vsitetop *vsitetop, int nvsitetop,
393 const char res[], const char atom1[],
394 const char atom2[], const char atom3[])
399 while (i < nvsitetop && gmx_strcasecmp(res, vsitetop[i].resname))
405 gmx_fatal(FARGS, "No vsite information for residue %s found in vsite database.\n", res);
408 while (j < vsitetop[i].nangles &&
409 ( strcmp(atom1, vsitetop[i].angle[j].atom1) ||
410 strcmp(atom2, vsitetop[i].angle[j].atom2) ||
411 strcmp(atom3, vsitetop[i].angle[j].atom3)) &&
412 ( strcmp(atom3, vsitetop[i].angle[j].atom1) ||
413 strcmp(atom2, vsitetop[i].angle[j].atom2) ||
414 strcmp(atom1, vsitetop[i].angle[j].atom3)))
418 if (j == vsitetop[i].nangles)
420 gmx_fatal(FARGS, "Couldnt find angle %s-%s-%s for residue %s in vsite database.\n", atom1, atom2, atom3, res);
423 return vsitetop[i].angle[j].value;
427 static void count_bonds(int atom, t_params *psb, char ***atomname,
428 int *nrbonds, int *nrHatoms, int Hatoms[], int *Heavy,
429 int *nrheavies, int heavies[])
431 int i, heavy, other, nrb, nrH, nrhv;
433 /* find heavy atom bound to this hydrogen */
435 for (i = 0; (i < psb->nr) && (heavy == NOTSET); i++)
437 if (psb->param[i].AI == atom)
439 heavy = psb->param[i].AJ;
441 else if (psb->param[i].AJ == atom)
443 heavy = psb->param[i].AI;
448 gmx_fatal(FARGS, "unbound hydrogen atom %d", atom+1);
450 /* find all atoms bound to heavy atom */
455 for (i = 0; i < psb->nr; i++)
457 if (psb->param[i].AI == heavy)
459 other = psb->param[i].AJ;
461 else if (psb->param[i].AJ == heavy)
463 other = psb->param[i].AI;
468 if (is_hydrogen(*(atomname[other])))
475 heavies[nrhv] = other;
487 static void print_bonds(FILE *fp, int o2n[],
488 int nrHatoms, int Hatoms[], int Heavy,
489 int nrheavies, int heavies[])
493 fprintf(fp, "Found: %d Hatoms: ", nrHatoms);
494 for (i = 0; i < nrHatoms; i++)
496 fprintf(fp, " %d", o2n[Hatoms[i]]+1);
498 fprintf(fp, "; %d Heavy atoms: %d", nrheavies+1, o2n[Heavy]+1);
499 for (i = 0; i < nrheavies; i++)
501 fprintf(fp, " %d", o2n[heavies[i]]+1);
506 static int get_atype(int atom, t_atoms *at, int nrtp, t_restp rtp[],
507 gmx_residuetype_t rt)
514 if (at->atom[atom].m)
516 type = at->atom[atom].type;
520 /* get type from rtp */
521 rtpp = get_restp(*(at->resinfo[at->atom[atom].resind].name), nrtp, rtp);
522 bNterm = gmx_residuetype_is_protein(rt, *(at->resinfo[at->atom[atom].resind].name)) &&
523 (at->atom[atom].resind == 0);
524 j = search_jtype(rtpp, *(at->atomname[atom]), bNterm);
525 type = rtpp->atom[j].type;
530 static int vsite_nm2type(const char *name, gpp_atomtype_t atype)
534 tp = get_atomtype_type(name, atype);
537 gmx_fatal(FARGS, "Dummy mass type (%s) not found in atom type database",
544 static real get_amass(int atom, t_atoms *at, int nrtp, t_restp rtp[],
545 gmx_residuetype_t rt)
552 if (at->atom[atom].m)
554 mass = at->atom[atom].m;
558 /* get mass from rtp */
559 rtpp = get_restp(*(at->resinfo[at->atom[atom].resind].name), nrtp, rtp);
560 bNterm = gmx_residuetype_is_protein(rt, *(at->resinfo[at->atom[atom].resind].name)) &&
561 (at->atom[atom].resind == 0);
562 j = search_jtype(rtpp, *(at->atomname[atom]), bNterm);
563 mass = rtpp->atom[j].m;
568 static void my_add_param(t_params *plist, int ai, int aj, real b)
570 static real c[MAXFORCEPARAM] =
571 { NOTSET, NOTSET, NOTSET, NOTSET, NOTSET, NOTSET };
574 add_param(plist, ai, aj, c, NULL);
577 static void add_vsites(t_params plist[], int vsite_type[],
578 int Heavy, int nrHatoms, int Hatoms[],
579 int nrheavies, int heavies[])
581 int i, j, ftype, other, moreheavy, bb;
582 gmx_bool bSwapParity;
584 for (i = 0; i < nrHatoms; i++)
586 ftype = vsite_type[Hatoms[i]];
587 /* Errors in setting the vsite_type should really be caugth earlier,
588 * because here it's not possible to print any useful error message.
589 * But it's still better to print a message than to segfault.
593 gmx_incons("Undetected error in setting up virtual sites");
595 bSwapParity = (ftype < 0);
596 vsite_type[Hatoms[i]] = ftype = abs(ftype);
597 if (ftype == F_BONDS)
599 if ( (nrheavies != 1) && (nrHatoms != 1) )
601 gmx_fatal(FARGS, "cannot make constraint in add_vsites for %d heavy "
602 "atoms and %d hydrogen atoms", nrheavies, nrHatoms);
604 my_add_param(&(plist[F_CONSTRNC]), Hatoms[i], heavies[0], NOTSET);
615 gmx_fatal(FARGS, "Not enough heavy atoms (%d) for %s (min 3)",
617 interaction_function[vsite_type[Hatoms[i]]].name);
619 add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], heavies[1],
626 moreheavy = heavies[1];
630 /* find more heavy atoms */
631 other = moreheavy = NOTSET;
632 for (j = 0; (j < plist[F_BONDS].nr) && (moreheavy == NOTSET); j++)
634 if (plist[F_BONDS].param[j].AI == heavies[0])
636 other = plist[F_BONDS].param[j].AJ;
638 else if (plist[F_BONDS].param[j].AJ == heavies[0])
640 other = plist[F_BONDS].param[j].AI;
642 if ( (other != NOTSET) && (other != Heavy) )
647 if (moreheavy == NOTSET)
649 gmx_fatal(FARGS, "Unbound molecule part %d-%d", Heavy+1, Hatoms[0]+1);
652 add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], moreheavy,
660 gmx_fatal(FARGS, "Not enough heavy atoms (%d) for %s (min 4)",
662 interaction_function[vsite_type[Hatoms[i]]].name);
664 add_vsite4_atoms(&plist[ftype],
665 Hatoms[0], Heavy, heavies[0], heavies[1], heavies[2]);
669 gmx_fatal(FARGS, "can't use add_vsites for interaction function %s",
670 interaction_function[vsite_type[Hatoms[i]]].name);
676 #define ANGLE_6RING (DEG2RAD*120)
678 /* cosine rule: a^2 = b^2 + c^2 - 2 b c cos(alpha) */
679 /* get a^2 when a, b and alpha are given: */
680 #define cosrule(b, c, alpha) ( sqr(b) + sqr(c) - 2*b*c*cos(alpha) )
681 /* get cos(alpha) when a, b and c are given: */
682 #define acosrule(a, b, c) ( (sqr(b)+sqr(c)-sqr(a))/(2*b*c) )
684 static int gen_vsites_6ring(t_atoms *at, int *vsite_type[], t_params plist[],
685 int nrfound, int *ats, real bond_cc, real bond_ch,
686 real xcom, gmx_bool bDoZ)
688 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
690 atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
695 real a, b, dCGCE, tmp1, tmp2, mtot, mG, mrest;
696 real xCG, yCG, xCE1, yCE1, xCE2, yCE2;
697 /* CG, CE1 and CE2 stay and each get a part of the total mass,
698 * so the c-o-m stays the same.
705 gmx_incons("Generating vsites on 6-rings");
709 /* constraints between CG, CE1 and CE2: */
710 dCGCE = sqrt( cosrule(bond_cc, bond_cc, ANGLE_6RING) );
711 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE);
712 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE2], dCGCE);
713 my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atCE2], dCGCE);
715 /* rest will be vsite3 */
718 for (i = 0; i < (bDoZ ? atNR : atHZ); i++)
720 mtot += at->atom[ats[i]].m;
721 if (i != atCG && i != atCE1 && i != atCE2 && (bDoZ || (i != atHZ && i != atCZ) ) )
723 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
724 (*vsite_type)[ats[i]] = F_VSITE3;
728 /* Distribute mass so center-of-mass stays the same.
729 * The center-of-mass in the call is defined with x=0 at
730 * the CE1-CE2 bond and y=0 at the line from CG to the middle of CE1-CE2 bond.
732 xCG = -bond_cc+bond_cc*cos(ANGLE_6RING);
735 yCE1 = bond_cc*sin(0.5*ANGLE_6RING);
737 yCE2 = -bond_cc*sin(0.5*ANGLE_6RING);
739 mG = at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = xcom*mtot/xCG;
741 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB =
742 at->atom[ats[atCE2]].m = at->atom[ats[atCE2]].mB = mrest / 2;
744 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
745 tmp1 = dCGCE*sin(ANGLE_6RING*0.5);
746 tmp2 = bond_cc*cos(0.5*ANGLE_6RING) + tmp1;
748 a = b = -bond_ch / tmp1;
750 add_vsite3_param(&plist[F_VSITE3],
751 ats[atHE1], ats[atCE1], ats[atCE2], ats[atCG], a, b);
752 add_vsite3_param(&plist[F_VSITE3],
753 ats[atHE2], ats[atCE2], ats[atCE1], ats[atCG], a, b);
754 /* CD1, CD2 and CZ: */
756 add_vsite3_param(&plist[F_VSITE3],
757 ats[atCD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
758 add_vsite3_param(&plist[F_VSITE3],
759 ats[atCD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
762 add_vsite3_param(&plist[F_VSITE3],
763 ats[atCZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
765 /* HD1, HD2 and HZ: */
766 a = b = ( bond_ch + tmp2 ) / tmp1;
767 add_vsite3_param(&plist[F_VSITE3],
768 ats[atHD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
769 add_vsite3_param(&plist[F_VSITE3],
770 ats[atHD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
773 add_vsite3_param(&plist[F_VSITE3],
774 ats[atHZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
780 static int gen_vsites_phe(t_atoms *at, int *vsite_type[], t_params plist[],
781 int nrfound, int *ats, t_vsitetop *vsitetop, int nvsitetop)
783 real bond_cc, bond_ch;
786 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
788 atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
791 real x[atNR], y[atNR];
792 /* Aromatic rings have 6-fold symmetry, so we only need one bond length.
793 * (angle is always 120 degrees).
795 bond_cc = get_ddb_bond(vsitetop, nvsitetop, "PHE", "CD1", "CE1");
796 bond_ch = get_ddb_bond(vsitetop, nvsitetop, "PHE", "CD1", "HD1");
798 x[atCG] = -bond_cc+bond_cc*cos(ANGLE_6RING);
801 y[atCD1] = bond_cc*sin(0.5*ANGLE_6RING);
802 x[atHD1] = x[atCD1]+bond_ch*cos(ANGLE_6RING);
803 y[atHD1] = y[atCD1]+bond_ch*sin(ANGLE_6RING);
806 x[atHE1] = x[atCE1]-bond_ch*cos(ANGLE_6RING);
807 y[atHE1] = y[atCE1]+bond_ch*sin(ANGLE_6RING);
809 y[atCD2] = -y[atCD1];
811 y[atHD2] = -y[atHD1];
813 y[atCE2] = -y[atCE1];
815 y[atHE2] = -y[atHE1];
816 x[atCZ] = bond_cc*cos(0.5*ANGLE_6RING);
818 x[atHZ] = x[atCZ]+bond_ch;
822 for (i = 0; i < atNR; i++)
824 xcom += x[i]*at->atom[ats[i]].m;
825 mtot += at->atom[ats[i]].m;
829 return gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, TRUE);
832 static void calc_vsite3_param(real xd, real yd, real xi, real yi, real xj, real yj,
833 real xk, real yk, real *a, real *b)
835 /* determine parameters by solving the equation system, since we know the
836 * virtual site coordinates here.
838 real dx_ij, dx_ik, dy_ij, dy_ik;
845 b_ij = sqrt(dx_ij*dx_ij+dy_ij*dy_ij);
846 b_ik = sqrt(dx_ik*dx_ik+dy_ik*dy_ik);
848 *a = ( (xd-xi)*dy_ik - dx_ik*(yd-yi) ) / (dx_ij*dy_ik - dx_ik*dy_ij);
849 *b = ( yd - yi - (*a)*dy_ij ) / dy_ik;
853 static int gen_vsites_trp(gpp_atomtype_t atype, rvec *newx[],
854 t_atom *newatom[], char ***newatomname[],
855 int *o2n[], int *newvsite_type[], int *newcgnr[],
856 t_symtab *symtab, int *nadd, rvec x[], int *cgnr[],
857 t_atoms *at, int *vsite_type[], t_params plist[],
858 int nrfound, int *ats, int add_shift,
859 t_vsitetop *vsitetop, int nvsitetop)
862 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
864 atCB, atCG, atCD1, atHD1, atCD2, atNE1, atHE1, atCE2, atCE3, atHE3,
865 atCZ2, atHZ2, atCZ3, atHZ3, atCH2, atHH2, atNR
867 /* weights for determining the COM's of both rings (M1 and M2): */
868 real mw[NMASS][atNR] = {
869 { 0, 1, 1, 1, 0.5, 1, 1, 0.5, 0, 0,
871 { 0, 0, 0, 0, 0.5, 0, 0, 0.5, 1, 1,
875 real xi[atNR], yi[atNR];
876 real xcom[NMASS], ycom[NMASS], I, alpha;
877 real lineA, lineB, dist;
878 real b_CD2_CE2, b_NE1_CE2, b_CG_CD2, b_CH2_HH2, b_CE2_CZ2;
879 real b_NE1_HE1, b_CD2_CE3, b_CE3_CZ3, b_CB_CG;
880 real b_CZ2_CH2, b_CZ2_HZ2, b_CD1_HD1, b_CE3_HE3;
881 real b_CG_CD1, b_CZ3_HZ3;
882 real a_NE1_CE2_CD2, a_CE2_CD2_CG, a_CB_CG_CD2, a_CE2_CD2_CE3;
883 real a_CB_CG_CD1, a_CD2_CG_CD1, a_CE2_CZ2_HZ2, a_CZ2_CH2_HH2;
884 real a_CD2_CE2_CZ2, a_CD2_CE3_CZ3, a_CE3_CZ3_HZ3, a_CG_CD1_HD1;
885 real a_CE2_CZ2_CH2, a_HE1_NE1_CE2, a_CD2_CE3_HE3;
887 int atM[NMASS], tpM, i, i0, j, nvsite;
888 real mwtot, mtot, mM[NMASS], dCBM1, dCBM2, dM1M2;
889 real a, b, c[MAXFORCEPARAM];
890 rvec r_ij, r_ik, t1, t2;
895 gmx_incons("atom types in gen_vsites_trp");
897 /* Get geometry from database */
898 b_CD2_CE2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CD2", "CE2");
899 b_NE1_CE2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "NE1", "CE2");
900 b_CG_CD1 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CG", "CD1");
901 b_CG_CD2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CG", "CD2");
902 b_CB_CG = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CB", "CG");
903 b_CE2_CZ2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CE2", "CZ2");
904 b_CD2_CE3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CD2", "CE3");
905 b_CE3_CZ3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CE3", "CZ3");
906 b_CZ2_CH2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CZ2", "CH2");
908 b_CD1_HD1 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CD1", "HD1");
909 b_CZ2_HZ2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CZ2", "HZ2");
910 b_NE1_HE1 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "NE1", "HE1");
911 b_CH2_HH2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CH2", "HH2");
912 b_CE3_HE3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CE3", "HE3");
913 b_CZ3_HZ3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CZ3", "HZ3");
915 a_NE1_CE2_CD2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "NE1", "CE2", "CD2");
916 a_CE2_CD2_CG = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CD2", "CG");
917 a_CB_CG_CD2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CB", "CG", "CD2");
918 a_CD2_CG_CD1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CG", "CD1");
919 a_CB_CG_CD1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CB", "CG", "CD1");
921 a_CE2_CD2_CE3 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CD2", "CE3");
922 a_CD2_CE2_CZ2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CE2", "CZ2");
923 a_CD2_CE3_CZ3 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CE3", "CZ3");
924 a_CE3_CZ3_HZ3 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE3", "CZ3", "HZ3");
925 a_CZ2_CH2_HH2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CZ2", "CH2", "HH2");
926 a_CE2_CZ2_HZ2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CZ2", "HZ2");
927 a_CE2_CZ2_CH2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CZ2", "CH2");
928 a_CG_CD1_HD1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CG", "CD1", "HD1");
929 a_HE1_NE1_CE2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "HE1", "NE1", "CE2");
930 a_CD2_CE3_HE3 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CE3", "HE3");
932 /* Calculate local coordinates.
933 * y-axis (x=0) is the bond CD2-CE2.
934 * x-axis (y=0) is perpendicular to the bond CD2-CE2 and
935 * intersects the middle of the bond.
938 yi[atCD2] = -0.5*b_CD2_CE2;
941 yi[atCE2] = 0.5*b_CD2_CE2;
943 xi[atNE1] = -b_NE1_CE2*sin(a_NE1_CE2_CD2);
944 yi[atNE1] = yi[atCE2]-b_NE1_CE2*cos(a_NE1_CE2_CD2);
946 xi[atCG] = -b_CG_CD2*sin(a_CE2_CD2_CG);
947 yi[atCG] = yi[atCD2]+b_CG_CD2*cos(a_CE2_CD2_CG);
949 alpha = a_CE2_CD2_CG + M_PI - a_CB_CG_CD2;
950 xi[atCB] = xi[atCG]-b_CB_CG*sin(alpha);
951 yi[atCB] = yi[atCG]+b_CB_CG*cos(alpha);
953 alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - M_PI;
954 xi[atCD1] = xi[atCG]-b_CG_CD1*sin(alpha);
955 yi[atCD1] = yi[atCG]+b_CG_CD1*cos(alpha);
957 xi[atCE3] = b_CD2_CE3*sin(a_CE2_CD2_CE3);
958 yi[atCE3] = yi[atCD2]+b_CD2_CE3*cos(a_CE2_CD2_CE3);
960 xi[atCZ2] = b_CE2_CZ2*sin(a_CD2_CE2_CZ2);
961 yi[atCZ2] = yi[atCE2]-b_CE2_CZ2*cos(a_CD2_CE2_CZ2);
963 alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - M_PI;
964 xi[atCZ3] = xi[atCE3]+b_CE3_CZ3*sin(alpha);
965 yi[atCZ3] = yi[atCE3]+b_CE3_CZ3*cos(alpha);
967 alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - M_PI;
968 xi[atCH2] = xi[atCZ2]+b_CZ2_CH2*sin(alpha);
969 yi[atCH2] = yi[atCZ2]-b_CZ2_CH2*cos(alpha);
972 alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - a_CG_CD1_HD1;
973 xi[atHD1] = xi[atCD1]-b_CD1_HD1*sin(alpha);
974 yi[atHD1] = yi[atCD1]+b_CD1_HD1*cos(alpha);
976 alpha = a_NE1_CE2_CD2 + M_PI - a_HE1_NE1_CE2;
977 xi[atHE1] = xi[atNE1]-b_NE1_HE1*sin(alpha);
978 yi[atHE1] = yi[atNE1]-b_NE1_HE1*cos(alpha);
980 alpha = a_CE2_CD2_CE3 + M_PI - a_CD2_CE3_HE3;
981 xi[atHE3] = xi[atCE3]+b_CE3_HE3*sin(alpha);
982 yi[atHE3] = yi[atCE3]+b_CE3_HE3*cos(alpha);
984 alpha = a_CD2_CE2_CZ2 + M_PI - a_CE2_CZ2_HZ2;
985 xi[atHZ2] = xi[atCZ2]+b_CZ2_HZ2*sin(alpha);
986 yi[atHZ2] = yi[atCZ2]-b_CZ2_HZ2*cos(alpha);
988 alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - a_CZ2_CH2_HH2;
989 xi[atHZ3] = xi[atCZ3]+b_CZ3_HZ3*sin(alpha);
990 yi[atHZ3] = yi[atCZ3]+b_CZ3_HZ3*cos(alpha);
992 alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - a_CE3_CZ3_HZ3;
993 xi[atHH2] = xi[atCH2]+b_CH2_HH2*sin(alpha);
994 yi[atHH2] = yi[atCH2]-b_CH2_HH2*cos(alpha);
996 /* Determine coeff. for the line CB-CG */
997 lineA = (yi[atCB]-yi[atCG])/(xi[atCB]-xi[atCG]);
998 lineB = yi[atCG]-lineA*xi[atCG];
1000 /* Calculate masses for each ring and put it on the dummy masses */
1001 for (j = 0; j < NMASS; j++)
1003 mM[j] = xcom[j] = ycom[j] = 0;
1005 for (i = 0; i < atNR; i++)
1009 for (j = 0; j < NMASS; j++)
1011 mM[j] += mw[j][i] * at->atom[ats[i]].m;
1012 xcom[j] += xi[i] * mw[j][i] * at->atom[ats[i]].m;
1013 ycom[j] += yi[i] * mw[j][i] * at->atom[ats[i]].m;
1017 for (j = 0; j < NMASS; j++)
1023 /* get dummy mass type */
1024 tpM = vsite_nm2type("MW", atype);
1025 /* make space for 2 masses: shift all atoms starting with CB */
1027 for (j = 0; j < NMASS; j++)
1029 atM[j] = i0+*nadd+j;
1033 fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, (*o2n)[i0]+1);
1036 for (j = i0; j < at->nr; j++)
1038 (*o2n)[j] = j+*nadd;
1040 srenew(*newx, at->nr+*nadd);
1041 srenew(*newatom, at->nr+*nadd);
1042 srenew(*newatomname, at->nr+*nadd);
1043 srenew(*newvsite_type, at->nr+*nadd);
1044 srenew(*newcgnr, at->nr+*nadd);
1045 for (j = 0; j < NMASS; j++)
1047 (*newatomname)[at->nr+*nadd-1-j] = NULL;
1050 /* Dummy masses will be placed at the center-of-mass in each ring. */
1052 /* calc initial position for dummy masses in real (non-local) coordinates.
1053 * Cheat by using the routine to calculate virtual site parameters. It is
1054 * much easier when we have the coordinates expressed in terms of
1057 rvec_sub(x[ats[atCB]], x[ats[atCG]], r_ij);
1058 rvec_sub(x[ats[atCD2]], x[ats[atCG]], r_ik);
1059 calc_vsite3_param(xcom[0], ycom[0], xi[atCG], yi[atCG], xi[atCB], yi[atCB],
1060 xi[atCD2], yi[atCD2], &a, &b);
1063 rvec_add(t1, t2, t1);
1064 rvec_add(t1, x[ats[atCG]], (*newx)[atM[0]]);
1066 calc_vsite3_param(xcom[1], ycom[1], xi[atCG], yi[atCG], xi[atCB], yi[atCB],
1067 xi[atCD2], yi[atCD2], &a, &b);
1070 rvec_add(t1, t2, t1);
1071 rvec_add(t1, x[ats[atCG]], (*newx)[atM[1]]);
1073 /* set parameters for the masses */
1074 for (j = 0; j < NMASS; j++)
1076 sprintf(name, "MW%d", j+1);
1077 (*newatomname) [atM[j]] = put_symtab(symtab, name);
1078 (*newatom) [atM[j]].m = (*newatom)[atM[j]].mB = mM[j];
1079 (*newatom) [atM[j]].q = (*newatom)[atM[j]].qB = 0.0;
1080 (*newatom) [atM[j]].type = (*newatom)[atM[j]].typeB = tpM;
1081 (*newatom) [atM[j]].ptype = eptAtom;
1082 (*newatom) [atM[j]].resind = at->atom[i0].resind;
1083 (*newvsite_type)[atM[j]] = NOTSET;
1084 (*newcgnr) [atM[j]] = (*cgnr)[i0];
1086 /* renumber cgnr: */
1087 for (i = i0; i < at->nr; i++)
1092 /* constraints between CB, M1 and M2 */
1093 /* 'add_shift' says which atoms won't be renumbered afterwards */
1094 dCBM1 = sqrt( sqr(xcom[0]-xi[atCB]) + sqr(ycom[0]-yi[atCB]) );
1095 dM1M2 = sqrt( sqr(xcom[0]-xcom[1]) + sqr(ycom[0]-ycom[1]) );
1096 dCBM2 = sqrt( sqr(xcom[1]-xi[atCB]) + sqr(ycom[1]-yi[atCB]) );
1097 my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift+atM[0], dCBM1);
1098 my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift+atM[1], dCBM2);
1099 my_add_param(&(plist[F_CONSTRNC]), add_shift+atM[0], add_shift+atM[1], dM1M2);
1101 /* rest will be vsite3 */
1103 for (i = 0; i < atNR; i++)
1107 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
1108 (*vsite_type)[ats[i]] = F_VSITE3;
1113 /* now define all vsites from M1, M2, CB, ie:
1114 r_d = r_M1 + a r_M1_M2 + b r_M1_CB */
1115 for (i = 0; i < atNR; i++)
1117 if ( (*vsite_type)[ats[i]] == F_VSITE3)
1119 calc_vsite3_param(xi[i], yi[i], xcom[0], ycom[0], xcom[1], ycom[1], xi[atCB], yi[atCB], &a, &b);
1120 add_vsite3_param(&plist[F_VSITE3],
1121 ats[i], add_shift+atM[0], add_shift+atM[1], ats[atCB], a, b);
1129 static int gen_vsites_tyr(gpp_atomtype_t atype, rvec *newx[],
1130 t_atom *newatom[], char ***newatomname[],
1131 int *o2n[], int *newvsite_type[], int *newcgnr[],
1132 t_symtab *symtab, int *nadd, rvec x[], int *cgnr[],
1133 t_atoms *at, int *vsite_type[], t_params plist[],
1134 int nrfound, int *ats, int add_shift,
1135 t_vsitetop *vsitetop, int nvsitetop)
1137 int nvsite, i, i0, j, atM, tpM;
1138 real dCGCE, dCEOH, dCGM, tmp1, a, b;
1139 real bond_cc, bond_ch, bond_co, bond_oh, angle_coh;
1141 real vmass, vdist, mM;
1145 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1147 atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
1148 atCZ, atOH, atHH, atNR
1150 real xi[atNR], yi[atNR];
1151 /* CG, CE1, CE2 (as in general 6-ring) and OH and HH stay,
1152 rest gets virtualized.
1153 Now we have two linked triangles with one improper keeping them flat */
1154 if (atNR != nrfound)
1156 gmx_incons("Number of atom types in gen_vsites_tyr");
1159 /* Aromatic rings have 6-fold symmetry, so we only need one bond length
1160 * for the ring part (angle is always 120 degrees).
1162 bond_cc = get_ddb_bond(vsitetop, nvsitetop, "TYR", "CD1", "CE1");
1163 bond_ch = get_ddb_bond(vsitetop, nvsitetop, "TYR", "CD1", "HD1");
1164 bond_co = get_ddb_bond(vsitetop, nvsitetop, "TYR", "CZ", "OH");
1165 bond_oh = get_ddb_bond(vsitetop, nvsitetop, "TYR", "OH", "HH");
1166 angle_coh = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TYR", "CZ", "OH", "HH");
1168 xi[atCG] = -bond_cc+bond_cc*cos(ANGLE_6RING);
1170 xi[atCD1] = -bond_cc;
1171 yi[atCD1] = bond_cc*sin(0.5*ANGLE_6RING);
1172 xi[atHD1] = xi[atCD1]+bond_ch*cos(ANGLE_6RING);
1173 yi[atHD1] = yi[atCD1]+bond_ch*sin(ANGLE_6RING);
1175 yi[atCE1] = yi[atCD1];
1176 xi[atHE1] = xi[atCE1]-bond_ch*cos(ANGLE_6RING);
1177 yi[atHE1] = yi[atCE1]+bond_ch*sin(ANGLE_6RING);
1178 xi[atCD2] = xi[atCD1];
1179 yi[atCD2] = -yi[atCD1];
1180 xi[atHD2] = xi[atHD1];
1181 yi[atHD2] = -yi[atHD1];
1182 xi[atCE2] = xi[atCE1];
1183 yi[atCE2] = -yi[atCE1];
1184 xi[atHE2] = xi[atHE1];
1185 yi[atHE2] = -yi[atHE1];
1186 xi[atCZ] = bond_cc*cos(0.5*ANGLE_6RING);
1188 xi[atOH] = xi[atCZ]+bond_co;
1192 for (i = 0; i < atOH; i++)
1194 xcom += xi[i]*at->atom[ats[i]].m;
1195 mtot += at->atom[ats[i]].m;
1199 /* first do 6 ring as default,
1200 except CZ (we'll do that different) and HZ (we don't have that): */
1201 nvsite = gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, FALSE);
1203 /* then construct CZ from the 2nd triangle */
1204 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
1205 a = b = 0.5 * bond_co / ( bond_co - bond_cc*cos(ANGLE_6RING) );
1206 add_vsite3_param(&plist[F_VSITE3],
1207 ats[atCZ], ats[atOH], ats[atCE1], ats[atCE2], a, b);
1208 at->atom[ats[atCZ]].m = at->atom[ats[atCZ]].mB = 0;
1210 /* constraints between CE1, CE2 and OH */
1211 dCGCE = sqrt( cosrule(bond_cc, bond_cc, ANGLE_6RING) );
1212 dCEOH = sqrt( cosrule(bond_cc, bond_co, ANGLE_6RING) );
1213 my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atOH], dCEOH);
1214 my_add_param(&(plist[F_CONSTRNC]), ats[atCE2], ats[atOH], dCEOH);
1216 /* We also want to constrain the angle C-O-H, but since CZ is constructed
1217 * we need to introduce a constraint to CG.
1218 * CG is much further away, so that will lead to instabilities in LINCS
1219 * when we constrain both CG-HH and OH-HH distances. Instead of requiring
1220 * the use of lincs_order=8 we introduce a dummy mass three times further
1221 * away from OH than HH. The mass is accordingly a third, with the remaining
1222 * 2/3 moved to OH. This shouldnt cause any problems since the forces will
1223 * apply to the HH constructed atom and not directly on the virtual mass.
1226 vdist = 2.0*bond_oh;
1227 mM = at->atom[ats[atHH]].m/2.0;
1228 at->atom[ats[atOH]].m += mM; /* add 1/2 of original H mass */
1229 at->atom[ats[atOH]].mB += mM; /* add 1/2 of original H mass */
1230 at->atom[ats[atHH]].m = at->atom[ats[atHH]].mB = 0;
1232 /* get dummy mass type */
1233 tpM = vsite_nm2type("MW", atype);
1234 /* make space for 1 mass: shift HH only */
1239 fprintf(stderr, "Inserting 1 dummy mass at %d\n", (*o2n)[i0]+1);
1242 for (j = i0; j < at->nr; j++)
1244 (*o2n)[j] = j+*nadd;
1246 srenew(*newx, at->nr+*nadd);
1247 srenew(*newatom, at->nr+*nadd);
1248 srenew(*newatomname, at->nr+*nadd);
1249 srenew(*newvsite_type, at->nr+*nadd);
1250 srenew(*newcgnr, at->nr+*nadd);
1251 (*newatomname)[at->nr+*nadd-1] = NULL;
1253 /* Calc the dummy mass initial position */
1254 rvec_sub(x[ats[atHH]], x[ats[atOH]], r1);
1256 rvec_add(r1, x[ats[atHH]], (*newx)[atM]);
1258 strcpy(name, "MW1");
1259 (*newatomname) [atM] = put_symtab(symtab, name);
1260 (*newatom) [atM].m = (*newatom)[atM].mB = mM;
1261 (*newatom) [atM].q = (*newatom)[atM].qB = 0.0;
1262 (*newatom) [atM].type = (*newatom)[atM].typeB = tpM;
1263 (*newatom) [atM].ptype = eptAtom;
1264 (*newatom) [atM].resind = at->atom[i0].resind;
1265 (*newvsite_type)[atM] = NOTSET;
1266 (*newcgnr) [atM] = (*cgnr)[i0];
1267 /* renumber cgnr: */
1268 for (i = i0; i < at->nr; i++)
1273 (*vsite_type)[ats[atHH]] = F_VSITE2;
1275 /* assume we also want the COH angle constrained: */
1276 tmp1 = bond_cc*cos(0.5*ANGLE_6RING) + dCGCE*sin(ANGLE_6RING*0.5) + bond_co;
1277 dCGM = sqrt( cosrule(tmp1, vdist, angle_coh) );
1278 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], add_shift+atM, dCGM);
1279 my_add_param(&(plist[F_CONSTRNC]), ats[atOH], add_shift+atM, vdist);
1281 add_vsite2_param(&plist[F_VSITE2],
1282 ats[atHH], ats[atOH], add_shift+atM, 1.0/2.0);
1286 static int gen_vsites_his(t_atoms *at, int *vsite_type[], t_params plist[],
1287 int nrfound, int *ats, t_vsitetop *vsitetop, int nvsitetop)
1290 real a, b, alpha, dCGCE1, dCGNE2;
1291 real sinalpha, cosalpha;
1292 real xcom, ycom, mtot;
1293 real mG, mrest, mCE1, mNE2;
1294 real b_CG_ND1, b_ND1_CE1, b_CE1_NE2, b_CG_CD2, b_CD2_NE2;
1295 real b_ND1_HD1, b_NE2_HE2, b_CE1_HE1, b_CD2_HD2;
1296 real a_CG_ND1_CE1, a_CG_CD2_NE2, a_ND1_CE1_NE2, a_CE1_NE2_CD2;
1297 real a_NE2_CE1_HE1, a_NE2_CD2_HD2, a_CE1_ND1_HD1, a_CE1_NE2_HE2;
1300 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1302 atCG, atND1, atHD1, atCD2, atHD2, atCE1, atHE1, atNE2, atHE2, atNR
1304 real x[atNR], y[atNR];
1306 /* CG, CE1 and NE2 stay, each gets part of the total mass,
1307 rest gets virtualized */
1308 /* check number of atoms, 3 hydrogens may be missing: */
1309 /* assert( nrfound >= atNR-3 || nrfound <= atNR );
1310 * Don't understand the above logic. Shouldn't it be && rather than || ???
1312 if ((nrfound < atNR-3) || (nrfound > atNR))
1314 gmx_incons("Generating vsites for HIS");
1317 /* avoid warnings about uninitialized variables */
1318 b_ND1_HD1 = b_NE2_HE2 = b_CE1_HE1 = b_CD2_HD2 = a_NE2_CE1_HE1 =
1319 a_NE2_CD2_HD2 = a_CE1_ND1_HD1 = a_CE1_NE2_HE2 = 0;
1321 if (ats[atHD1] != NOTSET)
1323 if (ats[atHE2] != NOTSET)
1325 sprintf(resname, "HISH");
1329 sprintf(resname, "HISA");
1334 sprintf(resname, "HISB");
1337 /* Get geometry from database */
1338 b_CG_ND1 = get_ddb_bond(vsitetop, nvsitetop, resname, "CG", "ND1");
1339 b_ND1_CE1 = get_ddb_bond(vsitetop, nvsitetop, resname, "ND1", "CE1");
1340 b_CE1_NE2 = get_ddb_bond(vsitetop, nvsitetop, resname, "CE1", "NE2");
1341 b_CG_CD2 = get_ddb_bond(vsitetop, nvsitetop, resname, "CG", "CD2");
1342 b_CD2_NE2 = get_ddb_bond(vsitetop, nvsitetop, resname, "CD2", "NE2");
1343 a_CG_ND1_CE1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CG", "ND1", "CE1");
1344 a_CG_CD2_NE2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CG", "CD2", "NE2");
1345 a_ND1_CE1_NE2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "ND1", "CE1", "NE2");
1346 a_CE1_NE2_CD2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CE1", "NE2", "CD2");
1348 if (ats[atHD1] != NOTSET)
1350 b_ND1_HD1 = get_ddb_bond(vsitetop, nvsitetop, resname, "ND1", "HD1");
1351 a_CE1_ND1_HD1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CE1", "ND1", "HD1");
1353 if (ats[atHE2] != NOTSET)
1355 b_NE2_HE2 = get_ddb_bond(vsitetop, nvsitetop, resname, "NE2", "HE2");
1356 a_CE1_NE2_HE2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CE1", "NE2", "HE2");
1358 if (ats[atHD2] != NOTSET)
1360 b_CD2_HD2 = get_ddb_bond(vsitetop, nvsitetop, resname, "CD2", "HD2");
1361 a_NE2_CD2_HD2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "NE2", "CD2", "HD2");
1363 if (ats[atHE1] != NOTSET)
1365 b_CE1_HE1 = get_ddb_bond(vsitetop, nvsitetop, resname, "CE1", "HE1");
1366 a_NE2_CE1_HE1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "NE2", "CE1", "HE1");
1369 /* constraints between CG, CE1 and NE1 */
1370 dCGCE1 = sqrt( cosrule(b_CG_ND1, b_ND1_CE1, a_CG_ND1_CE1) );
1371 dCGNE2 = sqrt( cosrule(b_CG_CD2, b_CD2_NE2, a_CG_CD2_NE2) );
1373 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE1);
1374 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atNE2], dCGNE2);
1375 /* we already have a constraint CE1-NE2, so we don't add it again */
1377 /* calculate the positions in a local frame of reference.
1378 * The x-axis is the line from CG that makes a right angle
1379 * with the bond CE1-NE2, and the y-axis the bond CE1-NE2.
1381 /* First calculate the x-axis intersection with y-axis (=yCE1).
1382 * Get cos(angle CG-CE1-NE2) :
1384 cosalpha = acosrule(dCGNE2, dCGCE1, b_CE1_NE2);
1386 y[atCE1] = cosalpha*dCGCE1;
1388 y[atNE2] = y[atCE1]-b_CE1_NE2;
1389 sinalpha = sqrt(1-cosalpha*cosalpha);
1390 x[atCG] = -sinalpha*dCGCE1;
1392 x[atHE1] = x[atHE2] = x[atHD1] = x[atHD2] = 0;
1393 y[atHE1] = y[atHE2] = y[atHD1] = y[atHD2] = 0;
1395 /* calculate ND1 and CD2 positions from CE1 and NE2 */
1397 x[atND1] = -b_ND1_CE1*sin(a_ND1_CE1_NE2);
1398 y[atND1] = y[atCE1]-b_ND1_CE1*cos(a_ND1_CE1_NE2);
1400 x[atCD2] = -b_CD2_NE2*sin(a_CE1_NE2_CD2);
1401 y[atCD2] = y[atNE2]+b_CD2_NE2*cos(a_CE1_NE2_CD2);
1403 /* And finally the hydrogen positions */
1404 if (ats[atHE1] != NOTSET)
1406 x[atHE1] = x[atCE1] + b_CE1_HE1*sin(a_NE2_CE1_HE1);
1407 y[atHE1] = y[atCE1] - b_CE1_HE1*cos(a_NE2_CE1_HE1);
1409 /* HD2 - first get (ccw) angle from (positive) y-axis */
1410 if (ats[atHD2] != NOTSET)
1412 alpha = a_CE1_NE2_CD2 + M_PI - a_NE2_CD2_HD2;
1413 x[atHD2] = x[atCD2] - b_CD2_HD2*sin(alpha);
1414 y[atHD2] = y[atCD2] + b_CD2_HD2*cos(alpha);
1416 if (ats[atHD1] != NOTSET)
1418 /* HD1 - first get (cw) angle from (positive) y-axis */
1419 alpha = a_ND1_CE1_NE2 + M_PI - a_CE1_ND1_HD1;
1420 x[atHD1] = x[atND1] - b_ND1_HD1*sin(alpha);
1421 y[atHD1] = y[atND1] - b_ND1_HD1*cos(alpha);
1423 if (ats[atHE2] != NOTSET)
1425 x[atHE2] = x[atNE2] + b_NE2_HE2*sin(a_CE1_NE2_HE2);
1426 y[atHE2] = y[atNE2] + b_NE2_HE2*cos(a_CE1_NE2_HE2);
1428 /* Have all coordinates now */
1430 /* calc center-of-mass; keep atoms CG, CE1, NE2 and
1431 * set the rest to vsite3
1433 mtot = xcom = ycom = 0;
1435 for (i = 0; i < atNR; i++)
1437 if (ats[i] != NOTSET)
1439 mtot += at->atom[ats[i]].m;
1440 xcom += x[i]*at->atom[ats[i]].m;
1441 ycom += y[i]*at->atom[ats[i]].m;
1442 if (i != atCG && i != atCE1 && i != atNE2)
1444 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
1445 (*vsite_type)[ats[i]] = F_VSITE3;
1450 if (nvsite+3 != nrfound)
1452 gmx_incons("Generating vsites for HIS");
1458 /* distribute mass so that com stays the same */
1459 mG = xcom*mtot/x[atCG];
1461 mCE1 = (ycom-y[atNE2])*mrest/(y[atCE1]-y[atNE2]);
1464 at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = mG;
1465 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB = mCE1;
1466 at->atom[ats[atNE2]].m = at->atom[ats[atNE2]].mB = mNE2;
1469 if (ats[atHE1] != NOTSET)
1471 calc_vsite3_param(x[atHE1], y[atHE1], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1472 x[atCG], y[atCG], &a, &b);
1473 add_vsite3_param(&plist[F_VSITE3],
1474 ats[atHE1], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1477 if (ats[atHE2] != NOTSET)
1479 calc_vsite3_param(x[atHE2], y[atHE2], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1480 x[atCG], y[atCG], &a, &b);
1481 add_vsite3_param(&plist[F_VSITE3],
1482 ats[atHE2], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1486 calc_vsite3_param(x[atND1], y[atND1], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1487 x[atCG], y[atCG], &a, &b);
1488 add_vsite3_param(&plist[F_VSITE3],
1489 ats[atND1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1492 calc_vsite3_param(x[atCD2], y[atCD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1493 x[atCG], y[atCG], &a, &b);
1494 add_vsite3_param(&plist[F_VSITE3],
1495 ats[atCD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1498 if (ats[atHD1] != NOTSET)
1500 calc_vsite3_param(x[atHD1], y[atHD1], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1501 x[atCG], y[atCG], &a, &b);
1502 add_vsite3_param(&plist[F_VSITE3],
1503 ats[atHD1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1506 if (ats[atHD2] != NOTSET)
1508 calc_vsite3_param(x[atHD2], y[atHD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1509 x[atCG], y[atCG], &a, &b);
1510 add_vsite3_param(&plist[F_VSITE3],
1511 ats[atHD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1516 static gmx_bool is_vsite(int vsite_type)
1518 if (vsite_type == NOTSET)
1522 switch (abs(vsite_type) )
1536 static char atomnamesuffix[] = "1234";
1538 void do_vsites(int nrtp, t_restp rtp[], gpp_atomtype_t atype,
1539 t_atoms *at, t_symtab *symtab, rvec *x[],
1540 t_params plist[], int *vsite_type[], int *cgnr[],
1541 real mHmult, gmx_bool bVsiteAromatics,
1544 #define MAXATOMSPERRESIDUE 16
1545 int i, j, k, m, i0, ni0, whatres, resind, add_shift, ftype, nvsite, nadd;
1547 int nrfound = 0, needed, nrbonds, nrHatoms, Heavy, nrheavies, tpM, tpHeavy;
1548 int Hatoms[4], heavies[4], bb;
1549 gmx_bool bWARNING, bAddVsiteParam, bFirstWater;
1551 gmx_bool *bResProcessed;
1552 real mHtot, mtot, fact, fact2;
1553 rvec rpar, rperp, temp;
1554 char name[10], tpname[32], nexttpname[32], *ch;
1556 int *o2n, *newvsite_type, *newcgnr, ats[MAXATOMSPERRESIDUE];
1559 char ***newatomname;
1563 int nvsiteconf, nvsitetop, cmplength;
1564 gmx_bool isN, planarN, bFound;
1565 gmx_residuetype_t rt;
1567 t_vsiteconf *vsiteconflist;
1568 /* pointer to a list of CH3/NH3/NH2 configuration entries.
1569 * See comments in read_vsite_database. It isnt beautiful,
1570 * but it had to be fixed, and I dont even want to try to
1571 * maintain this part of the code...
1573 t_vsitetop *vsitetop;
1574 /* Pointer to a list of geometry (bond/angle) entries for
1575 * residues like PHE, TRP, TYR, HIS, etc., where we need
1576 * to know the geometry to construct vsite aromatics.
1577 * Note that equilibrium geometry isnt necessarily the same
1578 * as the individual bond and angle values given in the
1579 * force field (rings can be strained).
1582 /* if bVsiteAromatics=TRUE do_vsites will specifically convert atoms in
1583 PHE, TRP, TYR and HIS to a construction of virtual sites */
1585 resPHE, resTRP, resTYR, resHIS, resNR
1587 const char *resnms[resNR] = { "PHE", "TRP", "TYR", "HIS" };
1588 /* Amber03 alternative names for termini */
1589 const char *resnmsN[resNR] = { "NPHE", "NTRP", "NTYR", "NHIS" };
1590 const char *resnmsC[resNR] = { "CPHE", "CTRP", "CTYR", "CHIS" };
1591 /* HIS can be known as HISH, HIS1, HISA, HID, HIE, HIP, etc. too */
1592 gmx_bool bPartial[resNR] = { FALSE, FALSE, FALSE, TRUE };
1593 /* the atnms for every residue MUST correspond to the enums in the
1594 gen_vsites_* (one for each residue) routines! */
1595 /* also the atom names in atnms MUST be in the same order as in the .rtp! */
1596 const char *atnms[resNR][MAXATOMSPERRESIDUE+1] = {
1598 "CD1", "HD1", "CD2", "HD2",
1599 "CE1", "HE1", "CE2", "HE2",
1603 "CD1", "HD1", "CD2",
1604 "NE1", "HE1", "CE2", "CE3", "HE3",
1605 "CZ2", "HZ2", "CZ3", "HZ3",
1606 "CH2", "HH2", NULL },
1608 "CD1", "HD1", "CD2", "HD2",
1609 "CE1", "HE1", "CE2", "HE2",
1610 "CZ", "OH", "HH", NULL },
1612 "ND1", "HD1", "CD2", "HD2",
1613 "CE1", "HE1", "NE2", "HE2", NULL }
1618 printf("Searching for atoms to make virtual sites ...\n");
1619 fprintf(debug, "# # # VSITES # # #\n");
1622 ndb = fflib_search_file_end(ffdir, ".vsd", FALSE, &db);
1624 vsiteconflist = NULL;
1627 for (f = 0; f < ndb; f++)
1629 read_vsite_database(db[f], &vsiteconflist, &nvsiteconf, &vsitetop, &nvsitetop);
1637 /* we need a marker for which atoms should *not* be renumbered afterwards */
1638 add_shift = 10*at->nr;
1639 /* make arrays where masses can be inserted into */
1641 snew(newatom, at->nr);
1642 snew(newatomname, at->nr);
1643 snew(newvsite_type, at->nr);
1644 snew(newcgnr, at->nr);
1645 /* make index array to tell where the atoms go to when masses are inserted */
1647 for (i = 0; i < at->nr; i++)
1651 /* make index to tell which residues were already processed */
1652 snew(bResProcessed, at->nres);
1654 gmx_residuetype_init(&rt);
1656 /* generate vsite constructions */
1657 /* loop over all atoms */
1659 for (i = 0; (i < at->nr); i++)
1661 if (at->atom[i].resind != resind)
1663 resind = at->atom[i].resind;
1664 resnm = *(at->resinfo[resind].name);
1666 /* first check for aromatics to virtualize */
1667 /* don't waste our effort on DNA, water etc. */
1668 /* Only do the vsite aromatic stuff when we reach the
1669 * CA atom, since there might be an X2/X3 group on the
1670 * N-terminus that must be treated first.
1672 if (bVsiteAromatics &&
1673 !strcmp(*(at->atomname[i]), "CA") &&
1674 !bResProcessed[resind] &&
1675 gmx_residuetype_is_protein(rt, *(at->resinfo[resind].name)) )
1677 /* mark this residue */
1678 bResProcessed[resind] = TRUE;
1679 /* find out if this residue needs converting */
1681 for (j = 0; j < resNR && whatres == NOTSET; j++)
1684 cmplength = bPartial[j] ? strlen(resnm)-1 : strlen(resnm);
1686 bFound = ((gmx_strncasecmp(resnm, resnms[j], cmplength) == 0) ||
1687 (gmx_strncasecmp(resnm, resnmsN[j], cmplength) == 0) ||
1688 (gmx_strncasecmp(resnm, resnmsC[j], cmplength) == 0));
1693 /* get atoms we will be needing for the conversion */
1695 for (k = 0; atnms[j][k]; k++)
1698 for (m = i; m < at->nr && at->atom[m].resind == resind && ats[k] == NOTSET; m++)
1700 if (gmx_strcasecmp(*(at->atomname[m]), atnms[j][k]) == 0)
1708 /* now k is number of atom names in atnms[j] */
1717 if (nrfound < needed)
1719 gmx_fatal(FARGS, "not enough atoms found (%d, need %d) in "
1720 "residue %s %d while\n "
1721 "generating aromatics virtual site construction",
1722 nrfound, needed, resnm, at->resinfo[resind].nr);
1724 /* Advance overall atom counter */
1728 /* the enums for every residue MUST correspond to atnms[residue] */
1734 fprintf(stderr, "PHE at %d\n", o2n[ats[0]]+1);
1736 nvsite += gen_vsites_phe(at, vsite_type, plist, nrfound, ats, vsitetop, nvsitetop);
1741 fprintf(stderr, "TRP at %d\n", o2n[ats[0]]+1);
1743 nvsite += gen_vsites_trp(atype, &newx, &newatom, &newatomname, &o2n,
1744 &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
1745 at, vsite_type, plist, nrfound, ats, add_shift, vsitetop, nvsitetop);
1750 fprintf(stderr, "TYR at %d\n", o2n[ats[0]]+1);
1752 nvsite += gen_vsites_tyr(atype, &newx, &newatom, &newatomname, &o2n,
1753 &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
1754 at, vsite_type, plist, nrfound, ats, add_shift, vsitetop, nvsitetop);
1759 fprintf(stderr, "HIS at %d\n", o2n[ats[0]]+1);
1761 nvsite += gen_vsites_his(at, vsite_type, plist, nrfound, ats, vsitetop, nvsitetop);
1764 /* this means this residue won't be processed */
1767 gmx_fatal(FARGS, "DEATH HORROR in do_vsites (%s:%d)",
1768 __FILE__, __LINE__);
1769 } /* switch whatres */
1770 /* skip back to beginning of residue */
1771 while (i > 0 && at->atom[i-1].resind == resind)
1775 } /* if bVsiteAromatics & is protein */
1777 /* now process the rest of the hydrogens */
1778 /* only process hydrogen atoms which are not already set */
1779 if ( ((*vsite_type)[i] == NOTSET) && is_hydrogen(*(at->atomname[i])))
1781 /* find heavy atom, count #bonds from it and #H atoms bound to it
1782 and return H atom numbers (Hatoms) and heavy atom numbers (heavies) */
1783 count_bonds(i, &plist[F_BONDS], at->atomname,
1784 &nrbonds, &nrHatoms, Hatoms, &Heavy, &nrheavies, heavies);
1785 /* get Heavy atom type */
1786 tpHeavy = get_atype(Heavy, at, nrtp, rtp, rt);
1787 strcpy(tpname, get_atomtype_name(tpHeavy, atype));
1790 bAddVsiteParam = TRUE;
1791 /* nested if's which check nrHatoms, nrbonds and atomname */
1797 (*vsite_type)[i] = F_BONDS;
1799 case 3: /* =CH-, -NH- or =NH+- */
1800 (*vsite_type)[i] = F_VSITE3FD;
1802 case 4: /* --CH- (tert) */
1803 /* The old type 4FD had stability issues, so
1804 * all new constructs should use 4FDN
1806 (*vsite_type)[i] = F_VSITE4FDN;
1808 /* Check parity of heavy atoms from coordinates */
1813 rvec_sub((*x)[aj], (*x)[ai], tmpmat[0]);
1814 rvec_sub((*x)[ak], (*x)[ai], tmpmat[1]);
1815 rvec_sub((*x)[al], (*x)[ai], tmpmat[2]);
1817 if (det(tmpmat) > 0)
1825 default: /* nrbonds != 2, 3 or 4 */
1830 else if ( /*(nrHatoms == 2) && (nrbonds == 2) && REMOVED this test
1832 (gmx_strncasecmp(*at->atomname[Heavy], "OW", 2) == 0) )
1834 bAddVsiteParam = FALSE; /* this is water: skip these hydrogens */
1837 bFirstWater = FALSE;
1841 "Not converting hydrogens in water to virtual sites\n");
1845 else if ( (nrHatoms == 2) && (nrbonds == 4) )
1847 /* -CH2- , -NH2+- */
1848 (*vsite_type)[Hatoms[0]] = F_VSITE3OUT;
1849 (*vsite_type)[Hatoms[1]] = -F_VSITE3OUT;
1853 /* 2 or 3 hydrogen atom, with 3 or 4 bonds in total to the heavy atom.
1854 * If it is a nitrogen, first check if it is planar.
1856 isN = planarN = FALSE;
1857 if ((nrHatoms == 2) && ((*at->atomname[Heavy])[0] == 'N'))
1860 j = nitrogen_is_planar(vsiteconflist, nvsiteconf, tpname);
1863 gmx_fatal(FARGS, "No vsite database NH2 entry for type %s\n", tpname);
1867 if ( (nrHatoms == 2) && (nrbonds == 3) && ( !isN || planarN ) )
1869 /* =CH2 or, if it is a nitrogen NH2, it is a planar one */
1870 (*vsite_type)[Hatoms[0]] = F_VSITE3FAD;
1871 (*vsite_type)[Hatoms[1]] = -F_VSITE3FAD;
1873 else if ( ( (nrHatoms == 2) && (nrbonds == 3) &&
1874 ( isN && !planarN ) ) ||
1875 ( (nrHatoms == 3) && (nrbonds == 4) ) )
1877 /* CH3, NH3 or non-planar NH2 group */
1878 int Hat_vsite_type[3] = { F_VSITE3, F_VSITE3OUT, F_VSITE3OUT };
1879 gmx_bool Hat_SwapParity[3] = { FALSE, TRUE, FALSE };
1883 fprintf(stderr, "-XH3 or nonplanar NH2 group at %d\n", i+1);
1885 bAddVsiteParam = FALSE; /* we'll do this ourselves! */
1886 /* -NH2 (umbrella), -NH3+ or -CH3 */
1887 (*vsite_type)[Heavy] = F_VSITE3;
1888 for (j = 0; j < nrHatoms; j++)
1890 (*vsite_type)[Hatoms[j]] = Hat_vsite_type[j];
1892 /* get dummy mass type from first char of heavy atom type (N or C) */
1894 strcpy(nexttpname, get_atomtype_name(get_atype(heavies[0], at, nrtp, rtp, rt), atype));
1895 ch = get_dummymass_name(vsiteconflist, nvsiteconf, tpname, nexttpname);
1901 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);
1905 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);
1913 tpM = vsite_nm2type(name, atype);
1914 /* make space for 2 masses: shift all atoms starting with 'Heavy' */
1920 fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, o2n[i0]+1);
1923 for (j = i0; j < at->nr; j++)
1928 srenew(newx, at->nr+nadd);
1929 srenew(newatom, at->nr+nadd);
1930 srenew(newatomname, at->nr+nadd);
1931 srenew(newvsite_type, at->nr+nadd);
1932 srenew(newcgnr, at->nr+nadd);
1934 for (j = 0; j < NMASS; j++)
1936 newatomname[at->nr+nadd-1-j] = NULL;
1939 /* calculate starting position for the masses */
1941 /* get atom masses, and set Heavy and Hatoms mass to zero */
1942 for (j = 0; j < nrHatoms; j++)
1944 mHtot += get_amass(Hatoms[j], at, nrtp, rtp, rt);
1945 at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
1947 mtot = mHtot + get_amass(Heavy, at, nrtp, rtp, rt);
1948 at->atom[Heavy].m = at->atom[Heavy].mB = 0;
1955 /* generate vectors parallel and perpendicular to rotational axis:
1956 * rpar = Heavy -> Hcom
1957 * rperp = Hcom -> H1 */
1959 for (j = 0; j < nrHatoms; j++)
1961 rvec_inc(rpar, (*x)[Hatoms[j]]);
1963 svmul(1.0/nrHatoms, rpar, rpar); /* rpar = ( H1+H2+H3 ) / 3 */
1964 rvec_dec(rpar, (*x)[Heavy]); /* - Heavy */
1965 rvec_sub((*x)[Hatoms[0]], (*x)[Heavy], rperp);
1966 rvec_dec(rperp, rpar); /* rperp = H1 - Heavy - rpar */
1967 /* calc mass positions */
1968 svmul(fact2, rpar, temp);
1969 for (j = 0; (j < NMASS); j++) /* xM = xN + fact2 * rpar +/- fact * rperp */
1971 rvec_add((*x)[Heavy], temp, newx[ni0+j]);
1973 svmul(fact, rperp, temp);
1974 rvec_inc(newx[ni0 ], temp);
1975 rvec_dec(newx[ni0+1], temp);
1976 /* set atom parameters for the masses */
1977 for (j = 0; (j < NMASS); j++)
1979 /* make name: "M??#" or "M?#" (? is atomname, # is number) */
1981 for (k = 0; (*at->atomname[Heavy])[k] && ( k < NMASS ); k++)
1983 name[k+1] = (*at->atomname[Heavy])[k];
1985 name[k+1] = atomnamesuffix[j];
1987 newatomname[ni0+j] = put_symtab(symtab, name);
1988 newatom[ni0+j].m = newatom[ni0+j].mB = mtot/NMASS;
1989 newatom[ni0+j].q = newatom[ni0+j].qB = 0.0;
1990 newatom[ni0+j].type = newatom[ni0+j].typeB = tpM;
1991 newatom[ni0+j].ptype = eptAtom;
1992 newatom[ni0+j].resind = at->atom[i0].resind;
1993 newvsite_type[ni0+j] = NOTSET;
1994 newcgnr[ni0+j] = (*cgnr)[i0];
1996 /* add constraints between dummy masses and to heavies[0] */
1997 /* 'add_shift' says which atoms won't be renumbered afterwards */
1998 my_add_param(&(plist[F_CONSTRNC]), heavies[0], add_shift+ni0, NOTSET);
1999 my_add_param(&(plist[F_CONSTRNC]), heavies[0], add_shift+ni0+1, NOTSET);
2000 my_add_param(&(plist[F_CONSTRNC]), add_shift+ni0, add_shift+ni0+1, NOTSET);
2002 /* generate Heavy, H1, H2 and H3 from M1, M2 and heavies[0] */
2003 /* note that vsite_type cannot be NOTSET, because we just set it */
2004 add_vsite3_atoms (&plist[(*vsite_type)[Heavy]],
2005 Heavy, heavies[0], add_shift+ni0, add_shift+ni0+1,
2007 for (j = 0; j < nrHatoms; j++)
2009 add_vsite3_atoms(&plist[(*vsite_type)[Hatoms[j]]],
2010 Hatoms[j], heavies[0], add_shift+ni0, add_shift+ni0+1,
2024 "Warning: cannot convert atom %d %s (bound to a heavy atom "
2026 " %d bonds and %d bound hydrogens atoms) to virtual site\n",
2027 i+1, *(at->atomname[i]), tpname, nrbonds, nrHatoms);
2031 /* add vsite parameters to topology,
2032 also get rid of negative vsite_types */
2033 add_vsites(plist, (*vsite_type), Heavy, nrHatoms, Hatoms,
2034 nrheavies, heavies);
2035 /* transfer mass of virtual site to Heavy atom */
2036 for (j = 0; j < nrHatoms; j++)
2038 if (is_vsite((*vsite_type)[Hatoms[j]]))
2040 at->atom[Heavy].m += at->atom[Hatoms[j]].m;
2041 at->atom[Heavy].mB = at->atom[Heavy].m;
2042 at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
2049 fprintf(debug, "atom %d: ", o2n[i]+1);
2050 print_bonds(debug, o2n, nrHatoms, Hatoms, Heavy, nrheavies, heavies);
2052 } /* if vsite NOTSET & is hydrogen */
2054 } /* for i < at->nr */
2056 gmx_residuetype_destroy(rt);
2060 fprintf(debug, "Before inserting new atoms:\n");
2061 for (i = 0; i < at->nr; i++)
2063 fprintf(debug, "%4d %4d %4s %4d %4s %6d %-10s\n", i+1, o2n[i]+1,
2064 at->atomname[i] ? *(at->atomname[i]) : "(NULL)",
2065 at->resinfo[at->atom[i].resind].nr,
2066 at->resinfo[at->atom[i].resind].name ?
2067 *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
2069 ((*vsite_type)[i] == NOTSET) ?
2070 "NOTSET" : interaction_function[(*vsite_type)[i]].name);
2072 fprintf(debug, "new atoms to be inserted:\n");
2073 for (i = 0; i < at->nr+nadd; i++)
2077 fprintf(debug, "%4d %4s %4d %6d %-10s\n", i+1,
2078 newatomname[i] ? *(newatomname[i]) : "(NULL)",
2079 newatom[i].resind, newcgnr[i],
2080 (newvsite_type[i] == NOTSET) ?
2081 "NOTSET" : interaction_function[newvsite_type[i]].name);
2086 /* add all original atoms to the new arrays, using o2n index array */
2087 for (i = 0; i < at->nr; i++)
2089 newatomname [o2n[i]] = at->atomname [i];
2090 newatom [o2n[i]] = at->atom [i];
2091 newvsite_type[o2n[i]] = (*vsite_type)[i];
2092 newcgnr [o2n[i]] = (*cgnr) [i];
2093 copy_rvec((*x)[i], newx[o2n[i]]);
2095 /* throw away old atoms */
2097 sfree(at->atomname);
2101 /* put in the new ones */
2104 at->atomname = newatomname;
2105 *vsite_type = newvsite_type;
2108 if (at->nr > add_shift)
2110 gmx_fatal(FARGS, "Added impossible amount of dummy masses "
2111 "(%d on a total of %d atoms)\n", nadd, at->nr-nadd);
2116 fprintf(debug, "After inserting new atoms:\n");
2117 for (i = 0; i < at->nr; i++)
2119 fprintf(debug, "%4d %4s %4d %4s %6d %-10s\n", i+1,
2120 at->atomname[i] ? *(at->atomname[i]) : "(NULL)",
2121 at->resinfo[at->atom[i].resind].nr,
2122 at->resinfo[at->atom[i].resind].name ?
2123 *(at->resinfo[at->atom[i].resind].name) : "(NULL)",
2125 ((*vsite_type)[i] == NOTSET) ?
2126 "NOTSET" : interaction_function[(*vsite_type)[i]].name);
2130 /* now renumber all the interactions because of the added atoms */
2131 for (ftype = 0; ftype < F_NRE; ftype++)
2133 params = &(plist[ftype]);
2136 fprintf(debug, "Renumbering %d %s\n", params->nr,
2137 interaction_function[ftype].longname);
2139 for (i = 0; i < params->nr; i++)
2141 for (j = 0; j < NRAL(ftype); j++)
2143 if (params->param[i].a[j] >= add_shift)
2147 fprintf(debug, " [%u -> %u]", params->param[i].a[j],
2148 params->param[i].a[j]-add_shift);
2150 params->param[i].a[j] = params->param[i].a[j]-add_shift;
2156 fprintf(debug, " [%u -> %d]", params->param[i].a[j],
2157 o2n[params->param[i].a[j]]);
2159 params->param[i].a[j] = o2n[params->param[i].a[j]];
2164 fprintf(debug, "\n");
2168 /* now check if atoms in the added constraints are in increasing order */
2169 params = &(plist[F_CONSTRNC]);
2170 for (i = 0; i < params->nr; i++)
2172 if (params->param[i].AI > params->param[i].AJ)
2174 j = params->param[i].AJ;
2175 params->param[i].AJ = params->param[i].AI;
2176 params->param[i].AI = j;
2183 /* tell the user what we did */
2184 fprintf(stderr, "Marked %d virtual sites\n", nvsite);
2185 fprintf(stderr, "Added %d dummy masses\n", nadd);
2186 fprintf(stderr, "Added %d new constraints\n", plist[F_CONSTRNC].nr);
2189 void do_h_mass(t_params *psb, int vsite_type[], t_atoms *at, real mHmult,
2190 gmx_bool bDeuterate)
2194 /* loop over all atoms */
2195 for (i = 0; i < at->nr; i++)
2197 /* adjust masses if i is hydrogen and not a virtual site */
2198 if (!is_vsite(vsite_type[i]) && is_hydrogen(*(at->atomname[i])) )
2200 /* find bonded heavy atom */
2202 for (j = 0; (j < psb->nr) && (a == NOTSET); j++)
2204 /* if other atom is not a virtual site, it is the one we want */
2205 if ( (psb->param[j].AI == i) &&
2206 !is_vsite(vsite_type[psb->param[j].AJ]) )
2208 a = psb->param[j].AJ;
2210 else if ( (psb->param[j].AJ == i) &&
2211 !is_vsite(vsite_type[psb->param[j].AI]) )
2213 a = psb->param[j].AI;
2218 gmx_fatal(FARGS, "Unbound hydrogen atom (%d) found while adjusting mass",
2222 /* adjust mass of i (hydrogen) with mHmult
2223 and correct mass of a (bonded atom) with same amount */
2226 at->atom[a].m -= (mHmult-1.0)*at->atom[i].m;
2227 at->atom[a].mB -= (mHmult-1.0)*at->atom[i].m;
2229 at->atom[i].m *= mHmult;
2230 at->atom[i].mB *= mHmult;