2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
39 #include "gen_vsite.h"
47 #include "gromacs/fileio/pdbio.h"
48 #include "gromacs/gmxpreprocess/add_par.h"
49 #include "gromacs/gmxpreprocess/fflibutil.h"
50 #include "gromacs/gmxpreprocess/gpp_atomtype.h"
51 #include "gromacs/gmxpreprocess/notset.h"
52 #include "gromacs/gmxpreprocess/resall.h"
53 #include "gromacs/gmxpreprocess/toputil.h"
54 #include "gromacs/math/functions.h"
55 #include "gromacs/math/units.h"
56 #include "gromacs/math/vec.h"
57 #include "gromacs/mdtypes/md_enums.h"
58 #include "gromacs/topology/ifunc.h"
59 #include "gromacs/topology/residuetypes.h"
60 #include "gromacs/topology/symtab.h"
61 #include "gromacs/utility/basedefinitions.h"
62 #include "gromacs/utility/cstringutil.h"
63 #include "gromacs/utility/fatalerror.h"
64 #include "gromacs/utility/futil.h"
65 #include "gromacs/utility/real.h"
66 #include "gromacs/utility/smalloc.h"
69 #define OPENDIR '[' /* starting sign for directive */
70 #define CLOSEDIR ']' /* ending sign for directive */
73 char atomtype
[MAXNAME
]; /* Type for the XH3/XH2 atom */
74 gmx_bool isplanar
; /* If true, the atomtype above and the three connected
75 * ones are in a planar geometry. The two next entries
76 * are undefined in that case
78 int nhydrogens
; /* number of connected hydrogens */
79 char nextheavytype
[MAXNAME
]; /* Type for the heavy atom bonded to XH2/XH3 */
80 char dummymass
[MAXNAME
]; /* The type of MNH* or MCH3* dummy mass to use */
84 /* Structure to represent average bond and angles values in vsite aromatic
85 * residues. Note that these are NOT necessarily the bonds and angles from the
86 * forcefield; many forcefields (like Amber, OPLS) have some inherent strain in
87 * 5-rings (i.e. the sum of angles is !=540, but impropers keep it planar)
90 char resname
[MAXNAME
];
93 struct vsitetop_bond
{
97 } *bond
; /* list of bonds */
98 struct vsitetop_angle
{
103 } *angle
; /* list of angles */
108 DDB_CH3
, DDB_NH3
, DDB_NH2
, DDB_PHE
, DDB_TYR
,
109 DDB_TRP
, DDB_HISA
, DDB_HISB
, DDB_HISH
, DDB_DIR_NR
112 typedef char t_dirname
[STRLEN
];
114 static const t_dirname ddb_dirnames
[DDB_DIR_NR
] = {
126 static int ddb_name2dir(char *name
)
128 /* Translate a directive name to the number of the directive.
129 * HID/HIE/HIP names are translated to the ones we use in Gromacs.
136 for (i
= 0; i
< DDB_DIR_NR
&& index
< 0; i
++)
138 if (!gmx_strcasecmp(name
, ddb_dirnames
[i
]))
148 static void read_vsite_database(const char *ddbname
,
149 t_vsiteconf
**pvsiteconflist
, int *nvsiteconf
,
150 t_vsitetop
**pvsitetoplist
, int *nvsitetop
)
152 /* This routine is a quick hack to fix the problem with hardcoded atomtypes
153 * and aromatic vsite parameters by reading them from a ff???.vsd file.
155 * The file can contain sections [ NH3 ], [ CH3 ], [ NH2 ], and ring residue names.
156 * For the NH3 and CH3 section each line has three fields. The first is the atomtype
157 * (nb: not bonded type) of the N/C atom to be replaced, the second field is
158 * the type of the next heavy atom it is bonded to, and the third field the type
159 * of dummy mass that will be used for this group.
161 * If the NH2 group planar (sp2 N) a different vsite construct is used, so in this
162 * case the second field should just be the word planar.
168 int i
, n
, k
, nvsite
, ntop
, curdir
;
169 t_vsiteconf
*vsiteconflist
;
170 t_vsitetop
*vsitetoplist
;
172 char s1
[MAXNAME
], s2
[MAXNAME
], s3
[MAXNAME
], s4
[MAXNAME
];
174 ddb
= libopen(ddbname
);
176 nvsite
= *nvsiteconf
;
177 vsiteconflist
= *pvsiteconflist
;
179 vsitetoplist
= *pvsitetoplist
;
183 snew(vsiteconflist
, 1);
184 snew(vsitetoplist
, 1);
186 while (fgets2(pline
, STRLEN
-2, ddb
) != nullptr)
188 strip_comment(pline
);
190 if (strlen(pline
) > 0)
192 if (pline
[0] == OPENDIR
)
194 strncpy(dirstr
, pline
+1, STRLEN
-2);
195 if ((ch
= strchr (dirstr
, CLOSEDIR
)) != nullptr)
201 if (!gmx_strcasecmp(dirstr
, "HID") ||
202 !gmx_strcasecmp(dirstr
, "HISD"))
204 sprintf(dirstr
, "HISA");
206 else if (!gmx_strcasecmp(dirstr
, "HIE") ||
207 !gmx_strcasecmp(dirstr
, "HISE"))
209 sprintf(dirstr
, "HISB");
211 else if (!gmx_strcasecmp(dirstr
, "HIP"))
213 sprintf(dirstr
, "HISH");
216 curdir
= ddb_name2dir(dirstr
);
219 gmx_fatal(FARGS
, "Invalid directive %s in vsite database %s",
228 gmx_fatal(FARGS
, "First entry in vsite database must be a directive.\n");
233 n
= sscanf(pline
, "%s%s%s", s1
, s2
, s3
);
234 if (n
< 3 && !gmx_strcasecmp(s2
, "planar"))
236 srenew(vsiteconflist
, nvsite
+1);
237 strncpy(vsiteconflist
[nvsite
].atomtype
, s1
, MAXNAME
-1);
238 vsiteconflist
[nvsite
].isplanar
= TRUE
;
239 vsiteconflist
[nvsite
].nextheavytype
[0] = 0;
240 vsiteconflist
[nvsite
].dummymass
[0] = 0;
241 vsiteconflist
[nvsite
].nhydrogens
= 2;
246 srenew(vsiteconflist
, (nvsite
+1));
247 strncpy(vsiteconflist
[nvsite
].atomtype
, s1
, MAXNAME
-1);
248 vsiteconflist
[nvsite
].isplanar
= FALSE
;
249 strncpy(vsiteconflist
[nvsite
].nextheavytype
, s2
, MAXNAME
-1);
250 strncpy(vsiteconflist
[nvsite
].dummymass
, s3
, MAXNAME
-1);
251 if (curdir
== DDB_NH2
)
253 vsiteconflist
[nvsite
].nhydrogens
= 2;
257 vsiteconflist
[nvsite
].nhydrogens
= 3;
263 gmx_fatal(FARGS
, "Not enough directives in vsite database line: %s\n", pline
);
273 while ((i
< ntop
) && gmx_strcasecmp(dirstr
, vsitetoplist
[i
].resname
))
277 /* Allocate a new topology entry if this is a new residue */
280 srenew(vsitetoplist
, ntop
+1);
281 ntop
++; /* i still points to current vsite topology entry */
282 strncpy(vsitetoplist
[i
].resname
, dirstr
, MAXNAME
-1);
283 vsitetoplist
[i
].nbonds
= vsitetoplist
[i
].nangles
= 0;
284 snew(vsitetoplist
[i
].bond
, 1);
285 snew(vsitetoplist
[i
].angle
, 1);
287 n
= sscanf(pline
, "%s%s%s%s", s1
, s2
, s3
, s4
);
291 k
= vsitetoplist
[i
].nbonds
++;
292 srenew(vsitetoplist
[i
].bond
, k
+1);
293 strncpy(vsitetoplist
[i
].bond
[k
].atom1
, s1
, MAXNAME
-1);
294 strncpy(vsitetoplist
[i
].bond
[k
].atom2
, s2
, MAXNAME
-1);
295 vsitetoplist
[i
].bond
[k
].value
= strtod(s3
, nullptr);
300 k
= vsitetoplist
[i
].nangles
++;
301 srenew(vsitetoplist
[i
].angle
, k
+1);
302 strncpy(vsitetoplist
[i
].angle
[k
].atom1
, s1
, MAXNAME
-1);
303 strncpy(vsitetoplist
[i
].angle
[k
].atom2
, s2
, MAXNAME
-1);
304 strncpy(vsitetoplist
[i
].angle
[k
].atom3
, s3
, MAXNAME
-1);
305 vsitetoplist
[i
].angle
[k
].value
= strtod(s4
, nullptr);
309 gmx_fatal(FARGS
, "Need 3 or 4 values to specify bond/angle values in %s: %s\n", ddbname
, pline
);
313 gmx_fatal(FARGS
, "Didnt find a case for directive %s in read_vsite_database\n", dirstr
);
319 *pvsiteconflist
= vsiteconflist
;
320 *pvsitetoplist
= vsitetoplist
;
321 *nvsiteconf
= nvsite
;
327 static int nitrogen_is_planar(t_vsiteconf vsiteconflist
[], int nvsiteconf
, char atomtype
[])
329 /* Return 1 if atomtype exists in database list and is planar, 0 if not,
330 * and -1 if not found.
333 gmx_bool found
= FALSE
;
334 for (i
= 0; i
< nvsiteconf
&& !found
; i
++)
336 found
= (!gmx_strcasecmp(vsiteconflist
[i
].atomtype
, atomtype
) && (vsiteconflist
[i
].nhydrogens
== 2));
340 res
= (vsiteconflist
[i
-1].isplanar
== TRUE
);
350 static char *get_dummymass_name(t_vsiteconf vsiteconflist
[], int nvsiteconf
, char atom
[], char nextheavy
[])
352 /* Return the dummy mass name if found, or NULL if not set in ddb database */
354 gmx_bool found
= FALSE
;
355 for (i
= 0; i
< nvsiteconf
&& !found
; i
++)
357 found
= (!gmx_strcasecmp(vsiteconflist
[i
].atomtype
, atom
) &&
358 !gmx_strcasecmp(vsiteconflist
[i
].nextheavytype
, nextheavy
));
362 return vsiteconflist
[i
-1].dummymass
;
372 static real
get_ddb_bond(t_vsitetop
*vsitetop
, int nvsitetop
,
374 const char atom1
[], const char atom2
[])
379 while (i
< nvsitetop
&& gmx_strcasecmp(res
, vsitetop
[i
].resname
))
385 gmx_fatal(FARGS
, "No vsite information for residue %s found in vsite database.\n", res
);
388 while (j
< vsitetop
[i
].nbonds
&&
389 ( strcmp(atom1
, vsitetop
[i
].bond
[j
].atom1
) || strcmp(atom2
, vsitetop
[i
].bond
[j
].atom2
)) &&
390 ( strcmp(atom2
, vsitetop
[i
].bond
[j
].atom1
) || strcmp(atom1
, vsitetop
[i
].bond
[j
].atom2
)))
394 if (j
== vsitetop
[i
].nbonds
)
396 gmx_fatal(FARGS
, "Couldnt find bond %s-%s for residue %s in vsite database.\n", atom1
, atom2
, res
);
399 return vsitetop
[i
].bond
[j
].value
;
403 static real
get_ddb_angle(t_vsitetop
*vsitetop
, int nvsitetop
,
404 const char res
[], const char atom1
[],
405 const char atom2
[], const char atom3
[])
410 while (i
< nvsitetop
&& gmx_strcasecmp(res
, vsitetop
[i
].resname
))
416 gmx_fatal(FARGS
, "No vsite information for residue %s found in vsite database.\n", res
);
419 while (j
< vsitetop
[i
].nangles
&&
420 ( strcmp(atom1
, vsitetop
[i
].angle
[j
].atom1
) ||
421 strcmp(atom2
, vsitetop
[i
].angle
[j
].atom2
) ||
422 strcmp(atom3
, vsitetop
[i
].angle
[j
].atom3
)) &&
423 ( strcmp(atom3
, vsitetop
[i
].angle
[j
].atom1
) ||
424 strcmp(atom2
, vsitetop
[i
].angle
[j
].atom2
) ||
425 strcmp(atom1
, vsitetop
[i
].angle
[j
].atom3
)))
429 if (j
== vsitetop
[i
].nangles
)
431 gmx_fatal(FARGS
, "Couldnt find angle %s-%s-%s for residue %s in vsite database.\n", atom1
, atom2
, atom3
, res
);
434 return vsitetop
[i
].angle
[j
].value
;
438 static void count_bonds(int atom
, t_params
*psb
, char ***atomname
,
439 int *nrbonds
, int *nrHatoms
, int Hatoms
[], int *Heavy
,
440 int *nrheavies
, int heavies
[])
442 int i
, heavy
, other
, nrb
, nrH
, nrhv
;
444 /* find heavy atom bound to this hydrogen */
446 for (i
= 0; (i
< psb
->nr
) && (heavy
== NOTSET
); i
++)
448 if (psb
->param
[i
].ai() == atom
)
450 heavy
= psb
->param
[i
].aj();
452 else if (psb
->param
[i
].aj() == atom
)
454 heavy
= psb
->param
[i
].ai();
459 gmx_fatal(FARGS
, "unbound hydrogen atom %d", atom
+1);
461 /* find all atoms bound to heavy atom */
466 for (i
= 0; i
< psb
->nr
; i
++)
468 if (psb
->param
[i
].ai() == heavy
)
470 other
= psb
->param
[i
].aj();
472 else if (psb
->param
[i
].aj() == heavy
)
474 other
= psb
->param
[i
].ai();
479 if (is_hydrogen(*(atomname
[other
])))
486 heavies
[nrhv
] = other
;
498 static void print_bonds(FILE *fp
, int o2n
[],
499 int nrHatoms
, int Hatoms
[], int Heavy
,
500 int nrheavies
, int heavies
[])
504 fprintf(fp
, "Found: %d Hatoms: ", nrHatoms
);
505 for (i
= 0; i
< nrHatoms
; i
++)
507 fprintf(fp
, " %d", o2n
[Hatoms
[i
]]+1);
509 fprintf(fp
, "; %d Heavy atoms: %d", nrheavies
+1, o2n
[Heavy
]+1);
510 for (i
= 0; i
< nrheavies
; i
++)
512 fprintf(fp
, " %d", o2n
[heavies
[i
]]+1);
517 static int get_atype(int atom
, t_atoms
*at
, int nrtp
, t_restp rtp
[],
518 gmx_residuetype_t
*rt
)
525 if (at
->atom
[atom
].m
)
527 type
= at
->atom
[atom
].type
;
531 /* get type from rtp */
532 rtpp
= get_restp(*(at
->resinfo
[at
->atom
[atom
].resind
].name
), nrtp
, rtp
);
533 bNterm
= gmx_residuetype_is_protein(rt
, *(at
->resinfo
[at
->atom
[atom
].resind
].name
)) &&
534 (at
->atom
[atom
].resind
== 0);
535 j
= search_jtype(rtpp
, *(at
->atomname
[atom
]), bNterm
);
536 type
= rtpp
->atom
[j
].type
;
541 static int vsite_nm2type(const char *name
, gpp_atomtype_t atype
)
545 tp
= get_atomtype_type(name
, atype
);
548 gmx_fatal(FARGS
, "Dummy mass type (%s) not found in atom type database",
555 static real
get_amass(int atom
, t_atoms
*at
, int nrtp
, t_restp rtp
[],
556 gmx_residuetype_t
*rt
)
563 if (at
->atom
[atom
].m
)
565 mass
= at
->atom
[atom
].m
;
569 /* get mass from rtp */
570 rtpp
= get_restp(*(at
->resinfo
[at
->atom
[atom
].resind
].name
), nrtp
, rtp
);
571 bNterm
= gmx_residuetype_is_protein(rt
, *(at
->resinfo
[at
->atom
[atom
].resind
].name
)) &&
572 (at
->atom
[atom
].resind
== 0);
573 j
= search_jtype(rtpp
, *(at
->atomname
[atom
]), bNterm
);
574 mass
= rtpp
->atom
[j
].m
;
579 static void my_add_param(t_params
*plist
, int ai
, int aj
, real b
)
581 static real c
[MAXFORCEPARAM
] =
582 { NOTSET
, NOTSET
, NOTSET
, NOTSET
, NOTSET
, NOTSET
};
585 add_param(plist
, ai
, aj
, c
, nullptr);
588 static void add_vsites(t_params plist
[], int vsite_type
[],
589 int Heavy
, int nrHatoms
, int Hatoms
[],
590 int nrheavies
, int heavies
[])
592 int i
, j
, ftype
, other
, moreheavy
;
593 gmx_bool bSwapParity
;
595 for (i
= 0; i
< nrHatoms
; i
++)
597 ftype
= vsite_type
[Hatoms
[i
]];
598 /* Errors in setting the vsite_type should really be caugth earlier,
599 * because here it's not possible to print any useful error message.
600 * But it's still better to print a message than to segfault.
604 gmx_incons("Undetected error in setting up virtual sites");
606 bSwapParity
= (ftype
< 0);
607 vsite_type
[Hatoms
[i
]] = ftype
= abs(ftype
);
608 if (ftype
== F_BONDS
)
610 if ( (nrheavies
!= 1) && (nrHatoms
!= 1) )
612 gmx_fatal(FARGS
, "cannot make constraint in add_vsites for %d heavy "
613 "atoms and %d hydrogen atoms", nrheavies
, nrHatoms
);
615 my_add_param(&(plist
[F_CONSTRNC
]), Hatoms
[i
], heavies
[0], NOTSET
);
626 gmx_fatal(FARGS
, "Not enough heavy atoms (%d) for %s (min 3)",
628 interaction_function
[vsite_type
[Hatoms
[i
]]].name
);
630 add_vsite3_atoms(&plist
[ftype
], Hatoms
[i
], Heavy
, heavies
[0], heavies
[1],
637 moreheavy
= heavies
[1];
641 /* find more heavy atoms */
642 other
= moreheavy
= NOTSET
;
643 for (j
= 0; (j
< plist
[F_BONDS
].nr
) && (moreheavy
== NOTSET
); j
++)
645 if (plist
[F_BONDS
].param
[j
].ai() == heavies
[0])
647 other
= plist
[F_BONDS
].param
[j
].aj();
649 else if (plist
[F_BONDS
].param
[j
].aj() == heavies
[0])
651 other
= plist
[F_BONDS
].param
[j
].ai();
653 if ( (other
!= NOTSET
) && (other
!= Heavy
) )
658 if (moreheavy
== NOTSET
)
660 gmx_fatal(FARGS
, "Unbound molecule part %d-%d", Heavy
+1, Hatoms
[0]+1);
663 add_vsite3_atoms(&plist
[ftype
], Hatoms
[i
], Heavy
, heavies
[0], moreheavy
,
671 gmx_fatal(FARGS
, "Not enough heavy atoms (%d) for %s (min 4)",
673 interaction_function
[vsite_type
[Hatoms
[i
]]].name
);
675 add_vsite4_atoms(&plist
[ftype
],
676 Hatoms
[0], Heavy
, heavies
[0], heavies
[1], heavies
[2]);
680 gmx_fatal(FARGS
, "can't use add_vsites for interaction function %s",
681 interaction_function
[vsite_type
[Hatoms
[i
]]].name
);
687 #define ANGLE_6RING (DEG2RAD*120)
689 /* cosine rule: a^2 = b^2 + c^2 - 2 b c cos(alpha) */
690 /* get a^2 when a, b and alpha are given: */
691 #define cosrule(b, c, alpha) ( gmx::square(b) + gmx::square(c) - 2*b*c*std::cos(alpha) )
692 /* get cos(alpha) when a, b and c are given: */
693 #define acosrule(a, b, c) ( (gmx::square(b)+gmx::square(c)-gmx::square(a))/(2*b*c) )
695 static int gen_vsites_6ring(t_atoms
*at
, int *vsite_type
[], t_params plist
[],
696 int nrfound
, int *ats
, real bond_cc
, real bond_ch
,
697 real xcom
, gmx_bool bDoZ
)
699 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
701 atCG
, atCD1
, atHD1
, atCD2
, atHD2
, atCE1
, atHE1
, atCE2
, atHE2
,
706 real a
, b
, dCGCE
, tmp1
, tmp2
, mtot
, mG
, mrest
;
708 /* CG, CE1 and CE2 stay and each get a part of the total mass,
709 * so the c-o-m stays the same.
716 gmx_incons("Generating vsites on 6-rings");
720 /* constraints between CG, CE1 and CE2: */
721 dCGCE
= std::sqrt( cosrule(bond_cc
, bond_cc
, ANGLE_6RING
) );
722 my_add_param(&(plist
[F_CONSTRNC
]), ats
[atCG
], ats
[atCE1
], dCGCE
);
723 my_add_param(&(plist
[F_CONSTRNC
]), ats
[atCG
], ats
[atCE2
], dCGCE
);
724 my_add_param(&(plist
[F_CONSTRNC
]), ats
[atCE1
], ats
[atCE2
], dCGCE
);
726 /* rest will be vsite3 */
729 for (i
= 0; i
< (bDoZ
? atNR
: atHZ
); i
++)
731 mtot
+= at
->atom
[ats
[i
]].m
;
732 if (i
!= atCG
&& i
!= atCE1
&& i
!= atCE2
&& (bDoZ
|| (i
!= atHZ
&& i
!= atCZ
) ) )
734 at
->atom
[ats
[i
]].m
= at
->atom
[ats
[i
]].mB
= 0;
735 (*vsite_type
)[ats
[i
]] = F_VSITE3
;
739 /* Distribute mass so center-of-mass stays the same.
740 * The center-of-mass in the call is defined with x=0 at
741 * the CE1-CE2 bond and y=0 at the line from CG to the middle of CE1-CE2 bond.
743 xCG
= -bond_cc
+bond_cc
*std::cos(ANGLE_6RING
);
745 mG
= at
->atom
[ats
[atCG
]].m
= at
->atom
[ats
[atCG
]].mB
= xcom
*mtot
/xCG
;
747 at
->atom
[ats
[atCE1
]].m
= at
->atom
[ats
[atCE1
]].mB
=
748 at
->atom
[ats
[atCE2
]].m
= at
->atom
[ats
[atCE2
]].mB
= mrest
/ 2;
750 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
751 tmp1
= dCGCE
*std::sin(ANGLE_6RING
*0.5);
752 tmp2
= bond_cc
*std::cos(0.5*ANGLE_6RING
) + tmp1
;
754 a
= b
= -bond_ch
/ tmp1
;
756 add_vsite3_param(&plist
[F_VSITE3
],
757 ats
[atHE1
], ats
[atCE1
], ats
[atCE2
], ats
[atCG
], a
, b
);
758 add_vsite3_param(&plist
[F_VSITE3
],
759 ats
[atHE2
], ats
[atCE2
], ats
[atCE1
], ats
[atCG
], a
, b
);
760 /* CD1, CD2 and CZ: */
762 add_vsite3_param(&plist
[F_VSITE3
],
763 ats
[atCD1
], ats
[atCE2
], ats
[atCE1
], ats
[atCG
], a
, b
);
764 add_vsite3_param(&plist
[F_VSITE3
],
765 ats
[atCD2
], ats
[atCE1
], ats
[atCE2
], ats
[atCG
], a
, b
);
768 add_vsite3_param(&plist
[F_VSITE3
],
769 ats
[atCZ
], ats
[atCG
], ats
[atCE1
], ats
[atCE2
], a
, b
);
771 /* HD1, HD2 and HZ: */
772 a
= b
= ( bond_ch
+ tmp2
) / tmp1
;
773 add_vsite3_param(&plist
[F_VSITE3
],
774 ats
[atHD1
], ats
[atCE2
], ats
[atCE1
], ats
[atCG
], a
, b
);
775 add_vsite3_param(&plist
[F_VSITE3
],
776 ats
[atHD2
], ats
[atCE1
], ats
[atCE2
], ats
[atCG
], a
, b
);
779 add_vsite3_param(&plist
[F_VSITE3
],
780 ats
[atHZ
], ats
[atCG
], ats
[atCE1
], ats
[atCE2
], a
, b
);
786 static int gen_vsites_phe(t_atoms
*at
, int *vsite_type
[], t_params plist
[],
787 int nrfound
, int *ats
, t_vsitetop
*vsitetop
, int nvsitetop
)
789 real bond_cc
, bond_ch
;
792 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
794 atCG
, atCD1
, atHD1
, atCD2
, atHD2
, atCE1
, atHE1
, atCE2
, atHE2
,
798 /* Aromatic rings have 6-fold symmetry, so we only need one bond length.
799 * (angle is always 120 degrees).
801 bond_cc
= get_ddb_bond(vsitetop
, nvsitetop
, "PHE", "CD1", "CE1");
802 bond_ch
= get_ddb_bond(vsitetop
, nvsitetop
, "PHE", "CD1", "HD1");
804 x
[atCG
] = -bond_cc
+bond_cc
*std::cos(ANGLE_6RING
);
806 x
[atHD1
] = x
[atCD1
]+bond_ch
*std::cos(ANGLE_6RING
);
808 x
[atHE1
] = x
[atCE1
]-bond_ch
*std::cos(ANGLE_6RING
);
813 x
[atCZ
] = bond_cc
*std::cos(0.5*ANGLE_6RING
);
814 x
[atHZ
] = x
[atCZ
]+bond_ch
;
817 for (i
= 0; i
< atNR
; i
++)
819 xcom
+= x
[i
]*at
->atom
[ats
[i
]].m
;
820 mtot
+= at
->atom
[ats
[i
]].m
;
824 return gen_vsites_6ring(at
, vsite_type
, plist
, nrfound
, ats
, bond_cc
, bond_ch
, xcom
, TRUE
);
827 static void calc_vsite3_param(real xd
, real yd
, real xi
, real yi
, real xj
, real yj
,
828 real xk
, real yk
, real
*a
, real
*b
)
830 /* determine parameters by solving the equation system, since we know the
831 * virtual site coordinates here.
833 real dx_ij
, dx_ik
, dy_ij
, dy_ik
;
840 *a
= ( (xd
-xi
)*dy_ik
- dx_ik
*(yd
-yi
) ) / (dx_ij
*dy_ik
- dx_ik
*dy_ij
);
841 *b
= ( yd
- yi
- (*a
)*dy_ij
) / dy_ik
;
845 static int gen_vsites_trp(gpp_atomtype_t atype
, rvec
*newx
[],
846 t_atom
*newatom
[], char ***newatomname
[],
847 int *o2n
[], int *newvsite_type
[], int *newcgnr
[],
848 t_symtab
*symtab
, int *nadd
, rvec x
[], int *cgnr
[],
849 t_atoms
*at
, int *vsite_type
[], t_params plist
[],
850 int nrfound
, int *ats
, int add_shift
,
851 t_vsitetop
*vsitetop
, int nvsitetop
)
854 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
856 atCB
, atCG
, atCD1
, atHD1
, atCD2
, atNE1
, atHE1
, atCE2
, atCE3
, atHE3
,
857 atCZ2
, atHZ2
, atCZ3
, atHZ3
, atCH2
, atHH2
, atNR
859 /* weights for determining the COM's of both rings (M1 and M2): */
860 real mw
[NMASS
][atNR
] = {
861 { 0, 1, 1, 1, 0.5, 1, 1, 0.5, 0, 0,
863 { 0, 0, 0, 0, 0.5, 0, 0, 0.5, 1, 1,
867 real xi
[atNR
], yi
[atNR
];
868 real xcom
[NMASS
], ycom
[NMASS
], alpha
;
869 real b_CD2_CE2
, b_NE1_CE2
, b_CG_CD2
, b_CH2_HH2
, b_CE2_CZ2
;
870 real b_NE1_HE1
, b_CD2_CE3
, b_CE3_CZ3
, b_CB_CG
;
871 real b_CZ2_CH2
, b_CZ2_HZ2
, b_CD1_HD1
, b_CE3_HE3
;
872 real b_CG_CD1
, b_CZ3_HZ3
;
873 real a_NE1_CE2_CD2
, a_CE2_CD2_CG
, a_CB_CG_CD2
, a_CE2_CD2_CE3
;
874 real a_CD2_CG_CD1
, a_CE2_CZ2_HZ2
, a_CZ2_CH2_HH2
;
875 real a_CD2_CE2_CZ2
, a_CD2_CE3_CZ3
, a_CE3_CZ3_HZ3
, a_CG_CD1_HD1
;
876 real a_CE2_CZ2_CH2
, a_HE1_NE1_CE2
, a_CD2_CE3_HE3
;
877 int atM
[NMASS
], tpM
, i
, i0
, j
, nvsite
;
878 real mM
[NMASS
], dCBM1
, dCBM2
, dM1M2
;
880 rvec r_ij
, r_ik
, t1
, t2
;
885 gmx_incons("atom types in gen_vsites_trp");
887 /* Get geometry from database */
888 b_CD2_CE2
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "CD2", "CE2");
889 b_NE1_CE2
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "NE1", "CE2");
890 b_CG_CD1
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "CG", "CD1");
891 b_CG_CD2
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "CG", "CD2");
892 b_CB_CG
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "CB", "CG");
893 b_CE2_CZ2
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "CE2", "CZ2");
894 b_CD2_CE3
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "CD2", "CE3");
895 b_CE3_CZ3
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "CE3", "CZ3");
896 b_CZ2_CH2
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "CZ2", "CH2");
898 b_CD1_HD1
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "CD1", "HD1");
899 b_CZ2_HZ2
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "CZ2", "HZ2");
900 b_NE1_HE1
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "NE1", "HE1");
901 b_CH2_HH2
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "CH2", "HH2");
902 b_CE3_HE3
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "CE3", "HE3");
903 b_CZ3_HZ3
= get_ddb_bond(vsitetop
, nvsitetop
, "TRP", "CZ3", "HZ3");
905 a_NE1_CE2_CD2
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "NE1", "CE2", "CD2");
906 a_CE2_CD2_CG
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "CE2", "CD2", "CG");
907 a_CB_CG_CD2
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "CB", "CG", "CD2");
908 a_CD2_CG_CD1
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "CD2", "CG", "CD1");
909 /*a_CB_CG_CD1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CB", "CG", "CD1"); unused */
911 a_CE2_CD2_CE3
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "CE2", "CD2", "CE3");
912 a_CD2_CE2_CZ2
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "CD2", "CE2", "CZ2");
913 a_CD2_CE3_CZ3
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "CD2", "CE3", "CZ3");
914 a_CE3_CZ3_HZ3
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "CE3", "CZ3", "HZ3");
915 a_CZ2_CH2_HH2
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "CZ2", "CH2", "HH2");
916 a_CE2_CZ2_HZ2
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "CE2", "CZ2", "HZ2");
917 a_CE2_CZ2_CH2
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "CE2", "CZ2", "CH2");
918 a_CG_CD1_HD1
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "CG", "CD1", "HD1");
919 a_HE1_NE1_CE2
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "HE1", "NE1", "CE2");
920 a_CD2_CE3_HE3
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TRP", "CD2", "CE3", "HE3");
922 /* Calculate local coordinates.
923 * y-axis (x=0) is the bond CD2-CE2.
924 * x-axis (y=0) is perpendicular to the bond CD2-CE2 and
925 * intersects the middle of the bond.
928 yi
[atCD2
] = -0.5*b_CD2_CE2
;
931 yi
[atCE2
] = 0.5*b_CD2_CE2
;
933 xi
[atNE1
] = -b_NE1_CE2
*std::sin(a_NE1_CE2_CD2
);
934 yi
[atNE1
] = yi
[atCE2
]-b_NE1_CE2
*std::cos(a_NE1_CE2_CD2
);
936 xi
[atCG
] = -b_CG_CD2
*std::sin(a_CE2_CD2_CG
);
937 yi
[atCG
] = yi
[atCD2
]+b_CG_CD2
*std::cos(a_CE2_CD2_CG
);
939 alpha
= a_CE2_CD2_CG
+ M_PI
- a_CB_CG_CD2
;
940 xi
[atCB
] = xi
[atCG
]-b_CB_CG
*std::sin(alpha
);
941 yi
[atCB
] = yi
[atCG
]+b_CB_CG
*std::cos(alpha
);
943 alpha
= a_CE2_CD2_CG
+ a_CD2_CG_CD1
- M_PI
;
944 xi
[atCD1
] = xi
[atCG
]-b_CG_CD1
*std::sin(alpha
);
945 yi
[atCD1
] = yi
[atCG
]+b_CG_CD1
*std::cos(alpha
);
947 xi
[atCE3
] = b_CD2_CE3
*std::sin(a_CE2_CD2_CE3
);
948 yi
[atCE3
] = yi
[atCD2
]+b_CD2_CE3
*std::cos(a_CE2_CD2_CE3
);
950 xi
[atCZ2
] = b_CE2_CZ2
*std::sin(a_CD2_CE2_CZ2
);
951 yi
[atCZ2
] = yi
[atCE2
]-b_CE2_CZ2
*std::cos(a_CD2_CE2_CZ2
);
953 alpha
= a_CE2_CD2_CE3
+ a_CD2_CE3_CZ3
- M_PI
;
954 xi
[atCZ3
] = xi
[atCE3
]+b_CE3_CZ3
*std::sin(alpha
);
955 yi
[atCZ3
] = yi
[atCE3
]+b_CE3_CZ3
*std::cos(alpha
);
957 alpha
= a_CD2_CE2_CZ2
+ a_CE2_CZ2_CH2
- M_PI
;
958 xi
[atCH2
] = xi
[atCZ2
]+b_CZ2_CH2
*std::sin(alpha
);
959 yi
[atCH2
] = yi
[atCZ2
]-b_CZ2_CH2
*std::cos(alpha
);
962 alpha
= a_CE2_CD2_CG
+ a_CD2_CG_CD1
- a_CG_CD1_HD1
;
963 xi
[atHD1
] = xi
[atCD1
]-b_CD1_HD1
*std::sin(alpha
);
964 yi
[atHD1
] = yi
[atCD1
]+b_CD1_HD1
*std::cos(alpha
);
966 alpha
= a_NE1_CE2_CD2
+ M_PI
- a_HE1_NE1_CE2
;
967 xi
[atHE1
] = xi
[atNE1
]-b_NE1_HE1
*std::sin(alpha
);
968 yi
[atHE1
] = yi
[atNE1
]-b_NE1_HE1
*std::cos(alpha
);
970 alpha
= a_CE2_CD2_CE3
+ M_PI
- a_CD2_CE3_HE3
;
971 xi
[atHE3
] = xi
[atCE3
]+b_CE3_HE3
*std::sin(alpha
);
972 yi
[atHE3
] = yi
[atCE3
]+b_CE3_HE3
*std::cos(alpha
);
974 alpha
= a_CD2_CE2_CZ2
+ M_PI
- a_CE2_CZ2_HZ2
;
975 xi
[atHZ2
] = xi
[atCZ2
]+b_CZ2_HZ2
*std::sin(alpha
);
976 yi
[atHZ2
] = yi
[atCZ2
]-b_CZ2_HZ2
*std::cos(alpha
);
978 alpha
= a_CD2_CE2_CZ2
+ a_CE2_CZ2_CH2
- a_CZ2_CH2_HH2
;
979 xi
[atHZ3
] = xi
[atCZ3
]+b_CZ3_HZ3
*std::sin(alpha
);
980 yi
[atHZ3
] = yi
[atCZ3
]+b_CZ3_HZ3
*std::cos(alpha
);
982 alpha
= a_CE2_CD2_CE3
+ a_CD2_CE3_CZ3
- a_CE3_CZ3_HZ3
;
983 xi
[atHH2
] = xi
[atCH2
]+b_CH2_HH2
*std::sin(alpha
);
984 yi
[atHH2
] = yi
[atCH2
]-b_CH2_HH2
*std::cos(alpha
);
986 /* Calculate masses for each ring and put it on the dummy masses */
987 for (j
= 0; j
< NMASS
; j
++)
989 mM
[j
] = xcom
[j
] = ycom
[j
] = 0;
991 for (i
= 0; i
< atNR
; i
++)
995 for (j
= 0; j
< NMASS
; j
++)
997 mM
[j
] += mw
[j
][i
] * at
->atom
[ats
[i
]].m
;
998 xcom
[j
] += xi
[i
] * mw
[j
][i
] * at
->atom
[ats
[i
]].m
;
999 ycom
[j
] += yi
[i
] * mw
[j
][i
] * at
->atom
[ats
[i
]].m
;
1003 for (j
= 0; j
< NMASS
; j
++)
1009 /* get dummy mass type */
1010 tpM
= vsite_nm2type("MW", atype
);
1011 /* make space for 2 masses: shift all atoms starting with CB */
1013 for (j
= 0; j
< NMASS
; j
++)
1015 atM
[j
] = i0
+*nadd
+j
;
1019 fprintf(stderr
, "Inserting %d dummy masses at %d\n", NMASS
, (*o2n
)[i0
]+1);
1022 for (j
= i0
; j
< at
->nr
; j
++)
1024 (*o2n
)[j
] = j
+*nadd
;
1026 srenew(*newx
, at
->nr
+*nadd
);
1027 srenew(*newatom
, at
->nr
+*nadd
);
1028 srenew(*newatomname
, at
->nr
+*nadd
);
1029 srenew(*newvsite_type
, at
->nr
+*nadd
);
1030 srenew(*newcgnr
, at
->nr
+*nadd
);
1031 for (j
= 0; j
< NMASS
; j
++)
1033 (*newatomname
)[at
->nr
+*nadd
-1-j
] = nullptr;
1036 /* Dummy masses will be placed at the center-of-mass in each ring. */
1038 /* calc initial position for dummy masses in real (non-local) coordinates.
1039 * Cheat by using the routine to calculate virtual site parameters. It is
1040 * much easier when we have the coordinates expressed in terms of
1043 rvec_sub(x
[ats
[atCB
]], x
[ats
[atCG
]], r_ij
);
1044 rvec_sub(x
[ats
[atCD2
]], x
[ats
[atCG
]], r_ik
);
1045 calc_vsite3_param(xcom
[0], ycom
[0], xi
[atCG
], yi
[atCG
], xi
[atCB
], yi
[atCB
],
1046 xi
[atCD2
], yi
[atCD2
], &a
, &b
);
1049 rvec_add(t1
, t2
, t1
);
1050 rvec_add(t1
, x
[ats
[atCG
]], (*newx
)[atM
[0]]);
1052 calc_vsite3_param(xcom
[1], ycom
[1], xi
[atCG
], yi
[atCG
], xi
[atCB
], yi
[atCB
],
1053 xi
[atCD2
], yi
[atCD2
], &a
, &b
);
1056 rvec_add(t1
, t2
, t1
);
1057 rvec_add(t1
, x
[ats
[atCG
]], (*newx
)[atM
[1]]);
1059 /* set parameters for the masses */
1060 for (j
= 0; j
< NMASS
; j
++)
1062 sprintf(name
, "MW%d", j
+1);
1063 (*newatomname
) [atM
[j
]] = put_symtab(symtab
, name
);
1064 (*newatom
) [atM
[j
]].m
= (*newatom
)[atM
[j
]].mB
= mM
[j
];
1065 (*newatom
) [atM
[j
]].q
= (*newatom
)[atM
[j
]].qB
= 0.0;
1066 (*newatom
) [atM
[j
]].type
= (*newatom
)[atM
[j
]].typeB
= tpM
;
1067 (*newatom
) [atM
[j
]].ptype
= eptAtom
;
1068 (*newatom
) [atM
[j
]].resind
= at
->atom
[i0
].resind
;
1069 (*newatom
) [atM
[j
]].elem
[0] = 'M';
1070 (*newatom
) [atM
[j
]].elem
[1] = '\0';
1071 (*newvsite_type
)[atM
[j
]] = NOTSET
;
1072 (*newcgnr
) [atM
[j
]] = (*cgnr
)[i0
];
1074 /* renumber cgnr: */
1075 for (i
= i0
; i
< at
->nr
; i
++)
1080 /* constraints between CB, M1 and M2 */
1081 /* 'add_shift' says which atoms won't be renumbered afterwards */
1082 dCBM1
= std::hypot( xcom
[0]-xi
[atCB
], ycom
[0]-yi
[atCB
] );
1083 dM1M2
= std::hypot( xcom
[0]-xcom
[1], ycom
[0]-ycom
[1] );
1084 dCBM2
= std::hypot( xcom
[1]-xi
[atCB
], ycom
[1]-yi
[atCB
] );
1085 my_add_param(&(plist
[F_CONSTRNC
]), ats
[atCB
], add_shift
+atM
[0], dCBM1
);
1086 my_add_param(&(plist
[F_CONSTRNC
]), ats
[atCB
], add_shift
+atM
[1], dCBM2
);
1087 my_add_param(&(plist
[F_CONSTRNC
]), add_shift
+atM
[0], add_shift
+atM
[1], dM1M2
);
1089 /* rest will be vsite3 */
1091 for (i
= 0; i
< atNR
; i
++)
1095 at
->atom
[ats
[i
]].m
= at
->atom
[ats
[i
]].mB
= 0;
1096 (*vsite_type
)[ats
[i
]] = F_VSITE3
;
1101 /* now define all vsites from M1, M2, CB, ie:
1102 r_d = r_M1 + a r_M1_M2 + b r_M1_CB */
1103 for (i
= 0; i
< atNR
; i
++)
1105 if ( (*vsite_type
)[ats
[i
]] == F_VSITE3
)
1107 calc_vsite3_param(xi
[i
], yi
[i
], xcom
[0], ycom
[0], xcom
[1], ycom
[1], xi
[atCB
], yi
[atCB
], &a
, &b
);
1108 add_vsite3_param(&plist
[F_VSITE3
],
1109 ats
[i
], add_shift
+atM
[0], add_shift
+atM
[1], ats
[atCB
], a
, b
);
1117 static int gen_vsites_tyr(gpp_atomtype_t atype
, rvec
*newx
[],
1118 t_atom
*newatom
[], char ***newatomname
[],
1119 int *o2n
[], int *newvsite_type
[], int *newcgnr
[],
1120 t_symtab
*symtab
, int *nadd
, rvec x
[], int *cgnr
[],
1121 t_atoms
*at
, int *vsite_type
[], t_params plist
[],
1122 int nrfound
, int *ats
, int add_shift
,
1123 t_vsitetop
*vsitetop
, int nvsitetop
)
1125 int nvsite
, i
, i0
, j
, atM
, tpM
;
1126 real dCGCE
, dCEOH
, dCGM
, tmp1
, a
, b
;
1127 real bond_cc
, bond_ch
, bond_co
, bond_oh
, angle_coh
;
1133 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1135 atCG
, atCD1
, atHD1
, atCD2
, atHD2
, atCE1
, atHE1
, atCE2
, atHE2
,
1136 atCZ
, atOH
, atHH
, atNR
1139 /* CG, CE1, CE2 (as in general 6-ring) and OH and HH stay,
1140 rest gets virtualized.
1141 Now we have two linked triangles with one improper keeping them flat */
1142 if (atNR
!= nrfound
)
1144 gmx_incons("Number of atom types in gen_vsites_tyr");
1147 /* Aromatic rings have 6-fold symmetry, so we only need one bond length
1148 * for the ring part (angle is always 120 degrees).
1150 bond_cc
= get_ddb_bond(vsitetop
, nvsitetop
, "TYR", "CD1", "CE1");
1151 bond_ch
= get_ddb_bond(vsitetop
, nvsitetop
, "TYR", "CD1", "HD1");
1152 bond_co
= get_ddb_bond(vsitetop
, nvsitetop
, "TYR", "CZ", "OH");
1153 bond_oh
= get_ddb_bond(vsitetop
, nvsitetop
, "TYR", "OH", "HH");
1154 angle_coh
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, "TYR", "CZ", "OH", "HH");
1156 xi
[atCG
] = -bond_cc
+bond_cc
*std::cos(ANGLE_6RING
);
1157 xi
[atCD1
] = -bond_cc
;
1158 xi
[atHD1
] = xi
[atCD1
]+bond_ch
*std::cos(ANGLE_6RING
);
1160 xi
[atHE1
] = xi
[atCE1
]-bond_ch
*std::cos(ANGLE_6RING
);
1161 xi
[atCD2
] = xi
[atCD1
];
1162 xi
[atHD2
] = xi
[atHD1
];
1163 xi
[atCE2
] = xi
[atCE1
];
1164 xi
[atHE2
] = xi
[atHE1
];
1165 xi
[atCZ
] = bond_cc
*std::cos(0.5*ANGLE_6RING
);
1166 xi
[atOH
] = xi
[atCZ
]+bond_co
;
1169 for (i
= 0; i
< atOH
; i
++)
1171 xcom
+= xi
[i
]*at
->atom
[ats
[i
]].m
;
1172 mtot
+= at
->atom
[ats
[i
]].m
;
1176 /* first do 6 ring as default,
1177 except CZ (we'll do that different) and HZ (we don't have that): */
1178 nvsite
= gen_vsites_6ring(at
, vsite_type
, plist
, nrfound
, ats
, bond_cc
, bond_ch
, xcom
, FALSE
);
1180 /* then construct CZ from the 2nd triangle */
1181 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
1182 a
= b
= 0.5 * bond_co
/ ( bond_co
- bond_cc
*std::cos(ANGLE_6RING
) );
1183 add_vsite3_param(&plist
[F_VSITE3
],
1184 ats
[atCZ
], ats
[atOH
], ats
[atCE1
], ats
[atCE2
], a
, b
);
1185 at
->atom
[ats
[atCZ
]].m
= at
->atom
[ats
[atCZ
]].mB
= 0;
1187 /* constraints between CE1, CE2 and OH */
1188 dCGCE
= std::sqrt( cosrule(bond_cc
, bond_cc
, ANGLE_6RING
) );
1189 dCEOH
= std::sqrt( cosrule(bond_cc
, bond_co
, ANGLE_6RING
) );
1190 my_add_param(&(plist
[F_CONSTRNC
]), ats
[atCE1
], ats
[atOH
], dCEOH
);
1191 my_add_param(&(plist
[F_CONSTRNC
]), ats
[atCE2
], ats
[atOH
], dCEOH
);
1193 /* We also want to constrain the angle C-O-H, but since CZ is constructed
1194 * we need to introduce a constraint to CG.
1195 * CG is much further away, so that will lead to instabilities in LINCS
1196 * when we constrain both CG-HH and OH-HH distances. Instead of requiring
1197 * the use of lincs_order=8 we introduce a dummy mass three times further
1198 * away from OH than HH. The mass is accordingly a third, with the remaining
1199 * 2/3 moved to OH. This shouldn't cause any problems since the forces will
1200 * apply to the HH constructed atom and not directly on the virtual mass.
1203 vdist
= 2.0*bond_oh
;
1204 mM
= at
->atom
[ats
[atHH
]].m
/2.0;
1205 at
->atom
[ats
[atOH
]].m
+= mM
; /* add 1/2 of original H mass */
1206 at
->atom
[ats
[atOH
]].mB
+= mM
; /* add 1/2 of original H mass */
1207 at
->atom
[ats
[atHH
]].m
= at
->atom
[ats
[atHH
]].mB
= 0;
1209 /* get dummy mass type */
1210 tpM
= vsite_nm2type("MW", atype
);
1211 /* make space for 1 mass: shift HH only */
1216 fprintf(stderr
, "Inserting 1 dummy mass at %d\n", (*o2n
)[i0
]+1);
1219 for (j
= i0
; j
< at
->nr
; j
++)
1221 (*o2n
)[j
] = j
+*nadd
;
1223 srenew(*newx
, at
->nr
+*nadd
);
1224 srenew(*newatom
, at
->nr
+*nadd
);
1225 srenew(*newatomname
, at
->nr
+*nadd
);
1226 srenew(*newvsite_type
, at
->nr
+*nadd
);
1227 srenew(*newcgnr
, at
->nr
+*nadd
);
1228 (*newatomname
)[at
->nr
+*nadd
-1] = nullptr;
1230 /* Calc the dummy mass initial position */
1231 rvec_sub(x
[ats
[atHH
]], x
[ats
[atOH
]], r1
);
1233 rvec_add(r1
, x
[ats
[atHH
]], (*newx
)[atM
]);
1235 strcpy(name
, "MW1");
1236 (*newatomname
) [atM
] = put_symtab(symtab
, name
);
1237 (*newatom
) [atM
].m
= (*newatom
)[atM
].mB
= mM
;
1238 (*newatom
) [atM
].q
= (*newatom
)[atM
].qB
= 0.0;
1239 (*newatom
) [atM
].type
= (*newatom
)[atM
].typeB
= tpM
;
1240 (*newatom
) [atM
].ptype
= eptAtom
;
1241 (*newatom
) [atM
].resind
= at
->atom
[i0
].resind
;
1242 (*newatom
) [atM
].elem
[0] = 'M';
1243 (*newatom
) [atM
].elem
[1] = '\0';
1244 (*newvsite_type
)[atM
] = NOTSET
;
1245 (*newcgnr
) [atM
] = (*cgnr
)[i0
];
1246 /* renumber cgnr: */
1247 for (i
= i0
; i
< at
->nr
; i
++)
1252 (*vsite_type
)[ats
[atHH
]] = F_VSITE2
;
1254 /* assume we also want the COH angle constrained: */
1255 tmp1
= bond_cc
*std::cos(0.5*ANGLE_6RING
) + dCGCE
*std::sin(ANGLE_6RING
*0.5) + bond_co
;
1256 dCGM
= std::sqrt( cosrule(tmp1
, vdist
, angle_coh
) );
1257 my_add_param(&(plist
[F_CONSTRNC
]), ats
[atCG
], add_shift
+atM
, dCGM
);
1258 my_add_param(&(plist
[F_CONSTRNC
]), ats
[atOH
], add_shift
+atM
, vdist
);
1260 add_vsite2_param(&plist
[F_VSITE2
],
1261 ats
[atHH
], ats
[atOH
], add_shift
+atM
, 1.0/2.0);
1265 static int gen_vsites_his(t_atoms
*at
, int *vsite_type
[], t_params plist
[],
1266 int nrfound
, int *ats
, t_vsitetop
*vsitetop
, int nvsitetop
)
1269 real a
, b
, alpha
, dCGCE1
, dCGNE2
;
1270 real sinalpha
, cosalpha
;
1271 real xcom
, ycom
, mtot
;
1272 real mG
, mrest
, mCE1
, mNE2
;
1273 real b_CG_ND1
, b_ND1_CE1
, b_CE1_NE2
, b_CG_CD2
, b_CD2_NE2
;
1274 real b_ND1_HD1
, b_NE2_HE2
, b_CE1_HE1
, b_CD2_HD2
;
1275 real a_CG_ND1_CE1
, a_CG_CD2_NE2
, a_ND1_CE1_NE2
, a_CE1_NE2_CD2
;
1276 real a_NE2_CE1_HE1
, a_NE2_CD2_HD2
, a_CE1_ND1_HD1
, a_CE1_NE2_HE2
;
1279 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1281 atCG
, atND1
, atHD1
, atCD2
, atHD2
, atCE1
, atHE1
, atNE2
, atHE2
, atNR
1283 real x
[atNR
], y
[atNR
];
1285 /* CG, CE1 and NE2 stay, each gets part of the total mass,
1286 rest gets virtualized */
1287 /* check number of atoms, 3 hydrogens may be missing: */
1288 /* assert( nrfound >= atNR-3 || nrfound <= atNR );
1289 * Don't understand the above logic. Shouldn't it be && rather than || ???
1291 if ((nrfound
< atNR
-3) || (nrfound
> atNR
))
1293 gmx_incons("Generating vsites for HIS");
1296 /* avoid warnings about uninitialized variables */
1297 b_ND1_HD1
= b_NE2_HE2
= b_CE1_HE1
= b_CD2_HD2
= a_NE2_CE1_HE1
=
1298 a_NE2_CD2_HD2
= a_CE1_ND1_HD1
= a_CE1_NE2_HE2
= 0;
1300 if (ats
[atHD1
] != NOTSET
)
1302 if (ats
[atHE2
] != NOTSET
)
1304 sprintf(resname
, "HISH");
1308 sprintf(resname
, "HISA");
1313 sprintf(resname
, "HISB");
1316 /* Get geometry from database */
1317 b_CG_ND1
= get_ddb_bond(vsitetop
, nvsitetop
, resname
, "CG", "ND1");
1318 b_ND1_CE1
= get_ddb_bond(vsitetop
, nvsitetop
, resname
, "ND1", "CE1");
1319 b_CE1_NE2
= get_ddb_bond(vsitetop
, nvsitetop
, resname
, "CE1", "NE2");
1320 b_CG_CD2
= get_ddb_bond(vsitetop
, nvsitetop
, resname
, "CG", "CD2");
1321 b_CD2_NE2
= get_ddb_bond(vsitetop
, nvsitetop
, resname
, "CD2", "NE2");
1322 a_CG_ND1_CE1
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, resname
, "CG", "ND1", "CE1");
1323 a_CG_CD2_NE2
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, resname
, "CG", "CD2", "NE2");
1324 a_ND1_CE1_NE2
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, resname
, "ND1", "CE1", "NE2");
1325 a_CE1_NE2_CD2
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, resname
, "CE1", "NE2", "CD2");
1327 if (ats
[atHD1
] != NOTSET
)
1329 b_ND1_HD1
= get_ddb_bond(vsitetop
, nvsitetop
, resname
, "ND1", "HD1");
1330 a_CE1_ND1_HD1
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, resname
, "CE1", "ND1", "HD1");
1332 if (ats
[atHE2
] != NOTSET
)
1334 b_NE2_HE2
= get_ddb_bond(vsitetop
, nvsitetop
, resname
, "NE2", "HE2");
1335 a_CE1_NE2_HE2
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, resname
, "CE1", "NE2", "HE2");
1337 if (ats
[atHD2
] != NOTSET
)
1339 b_CD2_HD2
= get_ddb_bond(vsitetop
, nvsitetop
, resname
, "CD2", "HD2");
1340 a_NE2_CD2_HD2
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, resname
, "NE2", "CD2", "HD2");
1342 if (ats
[atHE1
] != NOTSET
)
1344 b_CE1_HE1
= get_ddb_bond(vsitetop
, nvsitetop
, resname
, "CE1", "HE1");
1345 a_NE2_CE1_HE1
= DEG2RAD
*get_ddb_angle(vsitetop
, nvsitetop
, resname
, "NE2", "CE1", "HE1");
1348 /* constraints between CG, CE1 and NE1 */
1349 dCGCE1
= std::sqrt( cosrule(b_CG_ND1
, b_ND1_CE1
, a_CG_ND1_CE1
) );
1350 dCGNE2
= std::sqrt( cosrule(b_CG_CD2
, b_CD2_NE2
, a_CG_CD2_NE2
) );
1352 my_add_param(&(plist
[F_CONSTRNC
]), ats
[atCG
], ats
[atCE1
], dCGCE1
);
1353 my_add_param(&(plist
[F_CONSTRNC
]), ats
[atCG
], ats
[atNE2
], dCGNE2
);
1354 /* we already have a constraint CE1-NE2, so we don't add it again */
1356 /* calculate the positions in a local frame of reference.
1357 * The x-axis is the line from CG that makes a right angle
1358 * with the bond CE1-NE2, and the y-axis the bond CE1-NE2.
1360 /* First calculate the x-axis intersection with y-axis (=yCE1).
1361 * Get cos(angle CG-CE1-NE2) :
1363 cosalpha
= acosrule(dCGNE2
, dCGCE1
, b_CE1_NE2
);
1365 y
[atCE1
] = cosalpha
*dCGCE1
;
1367 y
[atNE2
] = y
[atCE1
]-b_CE1_NE2
;
1368 sinalpha
= std::sqrt(1-cosalpha
*cosalpha
);
1369 x
[atCG
] = -sinalpha
*dCGCE1
;
1371 x
[atHE1
] = x
[atHE2
] = x
[atHD1
] = x
[atHD2
] = 0;
1372 y
[atHE1
] = y
[atHE2
] = y
[atHD1
] = y
[atHD2
] = 0;
1374 /* calculate ND1 and CD2 positions from CE1 and NE2 */
1376 x
[atND1
] = -b_ND1_CE1
*std::sin(a_ND1_CE1_NE2
);
1377 y
[atND1
] = y
[atCE1
]-b_ND1_CE1
*std::cos(a_ND1_CE1_NE2
);
1379 x
[atCD2
] = -b_CD2_NE2
*std::sin(a_CE1_NE2_CD2
);
1380 y
[atCD2
] = y
[atNE2
]+b_CD2_NE2
*std::cos(a_CE1_NE2_CD2
);
1382 /* And finally the hydrogen positions */
1383 if (ats
[atHE1
] != NOTSET
)
1385 x
[atHE1
] = x
[atCE1
] + b_CE1_HE1
*std::sin(a_NE2_CE1_HE1
);
1386 y
[atHE1
] = y
[atCE1
] - b_CE1_HE1
*std::cos(a_NE2_CE1_HE1
);
1388 /* HD2 - first get (ccw) angle from (positive) y-axis */
1389 if (ats
[atHD2
] != NOTSET
)
1391 alpha
= a_CE1_NE2_CD2
+ M_PI
- a_NE2_CD2_HD2
;
1392 x
[atHD2
] = x
[atCD2
] - b_CD2_HD2
*std::sin(alpha
);
1393 y
[atHD2
] = y
[atCD2
] + b_CD2_HD2
*std::cos(alpha
);
1395 if (ats
[atHD1
] != NOTSET
)
1397 /* HD1 - first get (cw) angle from (positive) y-axis */
1398 alpha
= a_ND1_CE1_NE2
+ M_PI
- a_CE1_ND1_HD1
;
1399 x
[atHD1
] = x
[atND1
] - b_ND1_HD1
*std::sin(alpha
);
1400 y
[atHD1
] = y
[atND1
] - b_ND1_HD1
*std::cos(alpha
);
1402 if (ats
[atHE2
] != NOTSET
)
1404 x
[atHE2
] = x
[atNE2
] + b_NE2_HE2
*std::sin(a_CE1_NE2_HE2
);
1405 y
[atHE2
] = y
[atNE2
] + b_NE2_HE2
*std::cos(a_CE1_NE2_HE2
);
1407 /* Have all coordinates now */
1409 /* calc center-of-mass; keep atoms CG, CE1, NE2 and
1410 * set the rest to vsite3
1412 mtot
= xcom
= ycom
= 0;
1414 for (i
= 0; i
< atNR
; i
++)
1416 if (ats
[i
] != NOTSET
)
1418 mtot
+= at
->atom
[ats
[i
]].m
;
1419 xcom
+= x
[i
]*at
->atom
[ats
[i
]].m
;
1420 ycom
+= y
[i
]*at
->atom
[ats
[i
]].m
;
1421 if (i
!= atCG
&& i
!= atCE1
&& i
!= atNE2
)
1423 at
->atom
[ats
[i
]].m
= at
->atom
[ats
[i
]].mB
= 0;
1424 (*vsite_type
)[ats
[i
]] = F_VSITE3
;
1429 if (nvsite
+3 != nrfound
)
1431 gmx_incons("Generating vsites for HIS");
1437 /* distribute mass so that com stays the same */
1438 mG
= xcom
*mtot
/x
[atCG
];
1440 mCE1
= (ycom
-y
[atNE2
])*mrest
/(y
[atCE1
]-y
[atNE2
]);
1443 at
->atom
[ats
[atCG
]].m
= at
->atom
[ats
[atCG
]].mB
= mG
;
1444 at
->atom
[ats
[atCE1
]].m
= at
->atom
[ats
[atCE1
]].mB
= mCE1
;
1445 at
->atom
[ats
[atNE2
]].m
= at
->atom
[ats
[atNE2
]].mB
= mNE2
;
1448 if (ats
[atHE1
] != NOTSET
)
1450 calc_vsite3_param(x
[atHE1
], y
[atHE1
], x
[atCE1
], y
[atCE1
], x
[atNE2
], y
[atNE2
],
1451 x
[atCG
], y
[atCG
], &a
, &b
);
1452 add_vsite3_param(&plist
[F_VSITE3
],
1453 ats
[atHE1
], ats
[atCE1
], ats
[atNE2
], ats
[atCG
], a
, b
);
1456 if (ats
[atHE2
] != NOTSET
)
1458 calc_vsite3_param(x
[atHE2
], y
[atHE2
], x
[atNE2
], y
[atNE2
], x
[atCE1
], y
[atCE1
],
1459 x
[atCG
], y
[atCG
], &a
, &b
);
1460 add_vsite3_param(&plist
[F_VSITE3
],
1461 ats
[atHE2
], ats
[atNE2
], ats
[atCE1
], ats
[atCG
], a
, b
);
1465 calc_vsite3_param(x
[atND1
], y
[atND1
], x
[atNE2
], y
[atNE2
], x
[atCE1
], y
[atCE1
],
1466 x
[atCG
], y
[atCG
], &a
, &b
);
1467 add_vsite3_param(&plist
[F_VSITE3
],
1468 ats
[atND1
], ats
[atNE2
], ats
[atCE1
], ats
[atCG
], a
, b
);
1471 calc_vsite3_param(x
[atCD2
], y
[atCD2
], x
[atCE1
], y
[atCE1
], x
[atNE2
], y
[atNE2
],
1472 x
[atCG
], y
[atCG
], &a
, &b
);
1473 add_vsite3_param(&plist
[F_VSITE3
],
1474 ats
[atCD2
], ats
[atCE1
], ats
[atNE2
], ats
[atCG
], a
, b
);
1477 if (ats
[atHD1
] != NOTSET
)
1479 calc_vsite3_param(x
[atHD1
], y
[atHD1
], x
[atNE2
], y
[atNE2
], x
[atCE1
], y
[atCE1
],
1480 x
[atCG
], y
[atCG
], &a
, &b
);
1481 add_vsite3_param(&plist
[F_VSITE3
],
1482 ats
[atHD1
], ats
[atNE2
], ats
[atCE1
], ats
[atCG
], a
, b
);
1485 if (ats
[atHD2
] != NOTSET
)
1487 calc_vsite3_param(x
[atHD2
], y
[atHD2
], x
[atCE1
], y
[atCE1
], x
[atNE2
], y
[atNE2
],
1488 x
[atCG
], y
[atCG
], &a
, &b
);
1489 add_vsite3_param(&plist
[F_VSITE3
],
1490 ats
[atHD2
], ats
[atCE1
], ats
[atNE2
], ats
[atCG
], a
, b
);
1495 static gmx_bool
is_vsite(int vsite_type
)
1497 if (vsite_type
== NOTSET
)
1501 switch (abs(vsite_type
) )
1515 static char atomnamesuffix
[] = "1234";
1517 void do_vsites(int nrtp
, t_restp rtp
[], gpp_atomtype_t atype
,
1518 t_atoms
*at
, t_symtab
*symtab
, rvec
*x
[],
1519 t_params plist
[], int *vsite_type
[], int *cgnr
[],
1520 real mHmult
, gmx_bool bVsiteAromatics
,
1523 #define MAXATOMSPERRESIDUE 16
1524 int i
, j
, k
, m
, i0
, ni0
, whatres
, resind
, add_shift
, ftype
, nvsite
, nadd
;
1526 int nrfound
= 0, needed
, nrbonds
, nrHatoms
, Heavy
, nrheavies
, tpM
, tpHeavy
;
1527 int Hatoms
[4], heavies
[4];
1528 gmx_bool bWARNING
, bAddVsiteParam
, bFirstWater
;
1530 gmx_bool
*bResProcessed
;
1531 real mHtot
, mtot
, fact
, fact2
;
1532 rvec rpar
, rperp
, temp
;
1533 char name
[10], tpname
[32], nexttpname
[32], *ch
;
1535 int *o2n
, *newvsite_type
, *newcgnr
, ats
[MAXATOMSPERRESIDUE
];
1538 char ***newatomname
;
1539 char *resnm
= nullptr;
1542 int nvsiteconf
, nvsitetop
, cmplength
;
1543 gmx_bool isN
, planarN
, bFound
;
1544 gmx_residuetype_t
*rt
;
1546 t_vsiteconf
*vsiteconflist
;
1547 /* pointer to a list of CH3/NH3/NH2 configuration entries.
1548 * See comments in read_vsite_database. It isnt beautiful,
1549 * but it had to be fixed, and I dont even want to try to
1550 * maintain this part of the code...
1552 t_vsitetop
*vsitetop
;
1553 /* Pointer to a list of geometry (bond/angle) entries for
1554 * residues like PHE, TRP, TYR, HIS, etc., where we need
1555 * to know the geometry to construct vsite aromatics.
1556 * Note that equilibrium geometry isnt necessarily the same
1557 * as the individual bond and angle values given in the
1558 * force field (rings can be strained).
1561 /* if bVsiteAromatics=TRUE do_vsites will specifically convert atoms in
1562 PHE, TRP, TYR and HIS to a construction of virtual sites */
1564 resPHE
, resTRP
, resTYR
, resHIS
, resNR
1566 const char *resnms
[resNR
] = { "PHE", "TRP", "TYR", "HIS" };
1567 /* Amber03 alternative names for termini */
1568 const char *resnmsN
[resNR
] = { "NPHE", "NTRP", "NTYR", "NHIS" };
1569 const char *resnmsC
[resNR
] = { "CPHE", "CTRP", "CTYR", "CHIS" };
1570 /* HIS can be known as HISH, HIS1, HISA, HID, HIE, HIP, etc. too */
1571 gmx_bool bPartial
[resNR
] = { FALSE
, FALSE
, FALSE
, TRUE
};
1572 /* the atnms for every residue MUST correspond to the enums in the
1573 gen_vsites_* (one for each residue) routines! */
1574 /* also the atom names in atnms MUST be in the same order as in the .rtp! */
1575 const char *atnms
[resNR
][MAXATOMSPERRESIDUE
+1] = {
1577 "CD1", "HD1", "CD2", "HD2",
1578 "CE1", "HE1", "CE2", "HE2",
1579 "CZ", "HZ", nullptr },
1582 "CD1", "HD1", "CD2",
1583 "NE1", "HE1", "CE2", "CE3", "HE3",
1584 "CZ2", "HZ2", "CZ3", "HZ3",
1585 "CH2", "HH2", nullptr },
1587 "CD1", "HD1", "CD2", "HD2",
1588 "CE1", "HE1", "CE2", "HE2",
1589 "CZ", "OH", "HH", nullptr },
1591 "ND1", "HD1", "CD2", "HD2",
1592 "CE1", "HE1", "NE2", "HE2", nullptr }
1597 printf("Searching for atoms to make virtual sites ...\n");
1598 fprintf(debug
, "# # # VSITES # # #\n");
1601 ndb
= fflib_search_file_end(ffdir
, ".vsd", FALSE
, &db
);
1603 vsiteconflist
= nullptr;
1606 for (f
= 0; f
< ndb
; f
++)
1608 read_vsite_database(db
[f
], &vsiteconflist
, &nvsiteconf
, &vsitetop
, &nvsitetop
);
1616 /* we need a marker for which atoms should *not* be renumbered afterwards */
1617 add_shift
= 10*at
->nr
;
1618 /* make arrays where masses can be inserted into */
1620 snew(newatom
, at
->nr
);
1621 snew(newatomname
, at
->nr
);
1622 snew(newvsite_type
, at
->nr
);
1623 snew(newcgnr
, at
->nr
);
1624 /* make index array to tell where the atoms go to when masses are inserted */
1626 for (i
= 0; i
< at
->nr
; i
++)
1630 /* make index to tell which residues were already processed */
1631 snew(bResProcessed
, at
->nres
);
1633 gmx_residuetype_init(&rt
);
1635 /* generate vsite constructions */
1636 /* loop over all atoms */
1638 for (i
= 0; (i
< at
->nr
); i
++)
1640 if (at
->atom
[i
].resind
!= resind
)
1642 resind
= at
->atom
[i
].resind
;
1643 resnm
= *(at
->resinfo
[resind
].name
);
1645 /* first check for aromatics to virtualize */
1646 /* don't waste our effort on DNA, water etc. */
1647 /* Only do the vsite aromatic stuff when we reach the
1648 * CA atom, since there might be an X2/X3 group on the
1649 * N-terminus that must be treated first.
1651 if (bVsiteAromatics
&&
1652 !strcmp(*(at
->atomname
[i
]), "CA") &&
1653 !bResProcessed
[resind
] &&
1654 gmx_residuetype_is_protein(rt
, *(at
->resinfo
[resind
].name
)) )
1656 /* mark this residue */
1657 bResProcessed
[resind
] = TRUE
;
1658 /* find out if this residue needs converting */
1660 for (j
= 0; j
< resNR
&& whatres
== NOTSET
; j
++)
1663 cmplength
= bPartial
[j
] ? strlen(resnm
)-1 : strlen(resnm
);
1665 bFound
= ((gmx_strncasecmp(resnm
, resnms
[j
], cmplength
) == 0) ||
1666 (gmx_strncasecmp(resnm
, resnmsN
[j
], cmplength
) == 0) ||
1667 (gmx_strncasecmp(resnm
, resnmsC
[j
], cmplength
) == 0));
1672 /* get atoms we will be needing for the conversion */
1674 for (k
= 0; atnms
[j
][k
]; k
++)
1677 for (m
= i
; m
< at
->nr
&& at
->atom
[m
].resind
== resind
&& ats
[k
] == NOTSET
; m
++)
1679 if (gmx_strcasecmp(*(at
->atomname
[m
]), atnms
[j
][k
]) == 0)
1687 /* now k is number of atom names in atnms[j] */
1696 if (nrfound
< needed
)
1698 gmx_fatal(FARGS
, "not enough atoms found (%d, need %d) in "
1699 "residue %s %d while\n "
1700 "generating aromatics virtual site construction",
1701 nrfound
, needed
, resnm
, at
->resinfo
[resind
].nr
);
1703 /* Advance overall atom counter */
1707 /* the enums for every residue MUST correspond to atnms[residue] */
1713 fprintf(stderr
, "PHE at %d\n", o2n
[ats
[0]]+1);
1715 nvsite
+= gen_vsites_phe(at
, vsite_type
, plist
, nrfound
, ats
, vsitetop
, nvsitetop
);
1720 fprintf(stderr
, "TRP at %d\n", o2n
[ats
[0]]+1);
1722 nvsite
+= gen_vsites_trp(atype
, &newx
, &newatom
, &newatomname
, &o2n
,
1723 &newvsite_type
, &newcgnr
, symtab
, &nadd
, *x
, cgnr
,
1724 at
, vsite_type
, plist
, nrfound
, ats
, add_shift
, vsitetop
, nvsitetop
);
1729 fprintf(stderr
, "TYR at %d\n", o2n
[ats
[0]]+1);
1731 nvsite
+= gen_vsites_tyr(atype
, &newx
, &newatom
, &newatomname
, &o2n
,
1732 &newvsite_type
, &newcgnr
, symtab
, &nadd
, *x
, cgnr
,
1733 at
, vsite_type
, plist
, nrfound
, ats
, add_shift
, vsitetop
, nvsitetop
);
1738 fprintf(stderr
, "HIS at %d\n", o2n
[ats
[0]]+1);
1740 nvsite
+= gen_vsites_his(at
, vsite_type
, plist
, nrfound
, ats
, vsitetop
, nvsitetop
);
1743 /* this means this residue won't be processed */
1746 gmx_fatal(FARGS
, "DEATH HORROR in do_vsites (%s:%d)",
1747 __FILE__
, __LINE__
);
1748 } /* switch whatres */
1749 /* skip back to beginning of residue */
1750 while (i
> 0 && at
->atom
[i
-1].resind
== resind
)
1754 } /* if bVsiteAromatics & is protein */
1756 /* now process the rest of the hydrogens */
1757 /* only process hydrogen atoms which are not already set */
1758 if ( ((*vsite_type
)[i
] == NOTSET
) && is_hydrogen(*(at
->atomname
[i
])))
1760 /* find heavy atom, count #bonds from it and #H atoms bound to it
1761 and return H atom numbers (Hatoms) and heavy atom numbers (heavies) */
1762 count_bonds(i
, &plist
[F_BONDS
], at
->atomname
,
1763 &nrbonds
, &nrHatoms
, Hatoms
, &Heavy
, &nrheavies
, heavies
);
1764 /* get Heavy atom type */
1765 tpHeavy
= get_atype(Heavy
, at
, nrtp
, rtp
, rt
);
1766 strcpy(tpname
, get_atomtype_name(tpHeavy
, atype
));
1769 bAddVsiteParam
= TRUE
;
1770 /* nested if's which check nrHatoms, nrbonds and atomname */
1776 (*vsite_type
)[i
] = F_BONDS
;
1778 case 3: /* =CH-, -NH- or =NH+- */
1779 (*vsite_type
)[i
] = F_VSITE3FD
;
1781 case 4: /* --CH- (tert) */
1782 /* The old type 4FD had stability issues, so
1783 * all new constructs should use 4FDN
1785 (*vsite_type
)[i
] = F_VSITE4FDN
;
1787 /* Check parity of heavy atoms from coordinates */
1792 rvec_sub((*x
)[aj
], (*x
)[ai
], tmpmat
[0]);
1793 rvec_sub((*x
)[ak
], (*x
)[ai
], tmpmat
[1]);
1794 rvec_sub((*x
)[al
], (*x
)[ai
], tmpmat
[2]);
1796 if (det(tmpmat
) > 0)
1804 default: /* nrbonds != 2, 3 or 4 */
1809 else if ( (nrHatoms
== 2) && (nrbonds
== 2) &&
1810 (at
->atom
[Heavy
].atomnumber
== 8) )
1812 bAddVsiteParam
= FALSE
; /* this is water: skip these hydrogens */
1815 bFirstWater
= FALSE
;
1819 "Not converting hydrogens in water to virtual sites\n");
1823 else if ( (nrHatoms
== 2) && (nrbonds
== 4) )
1825 /* -CH2- , -NH2+- */
1826 (*vsite_type
)[Hatoms
[0]] = F_VSITE3OUT
;
1827 (*vsite_type
)[Hatoms
[1]] = -F_VSITE3OUT
;
1831 /* 2 or 3 hydrogen atom, with 3 or 4 bonds in total to the heavy atom.
1832 * If it is a nitrogen, first check if it is planar.
1834 isN
= planarN
= FALSE
;
1835 if ((nrHatoms
== 2) && ((*at
->atomname
[Heavy
])[0] == 'N'))
1838 j
= nitrogen_is_planar(vsiteconflist
, nvsiteconf
, tpname
);
1841 gmx_fatal(FARGS
, "No vsite database NH2 entry for type %s\n", tpname
);
1845 if ( (nrHatoms
== 2) && (nrbonds
== 3) && ( !isN
|| planarN
) )
1847 /* =CH2 or, if it is a nitrogen NH2, it is a planar one */
1848 (*vsite_type
)[Hatoms
[0]] = F_VSITE3FAD
;
1849 (*vsite_type
)[Hatoms
[1]] = -F_VSITE3FAD
;
1851 else if ( ( (nrHatoms
== 2) && (nrbonds
== 3) &&
1852 ( isN
&& !planarN
) ) ||
1853 ( (nrHatoms
== 3) && (nrbonds
== 4) ) )
1855 /* CH3, NH3 or non-planar NH2 group */
1856 int Hat_vsite_type
[3] = { F_VSITE3
, F_VSITE3OUT
, F_VSITE3OUT
};
1857 gmx_bool Hat_SwapParity
[3] = { FALSE
, TRUE
, FALSE
};
1861 fprintf(stderr
, "-XH3 or nonplanar NH2 group at %d\n", i
+1);
1863 bAddVsiteParam
= FALSE
; /* we'll do this ourselves! */
1864 /* -NH2 (umbrella), -NH3+ or -CH3 */
1865 (*vsite_type
)[Heavy
] = F_VSITE3
;
1866 for (j
= 0; j
< nrHatoms
; j
++)
1868 (*vsite_type
)[Hatoms
[j
]] = Hat_vsite_type
[j
];
1870 /* get dummy mass type from first char of heavy atom type (N or C) */
1872 strcpy(nexttpname
, get_atomtype_name(get_atype(heavies
[0], at
, nrtp
, rtp
, rt
), atype
));
1873 ch
= get_dummymass_name(vsiteconflist
, nvsiteconf
, tpname
, nexttpname
);
1879 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
);
1883 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
);
1891 tpM
= vsite_nm2type(name
, atype
);
1892 /* make space for 2 masses: shift all atoms starting with 'Heavy' */
1898 fprintf(stderr
, "Inserting %d dummy masses at %d\n", NMASS
, o2n
[i0
]+1);
1901 for (j
= i0
; j
< at
->nr
; j
++)
1906 srenew(newx
, at
->nr
+nadd
);
1907 srenew(newatom
, at
->nr
+nadd
);
1908 srenew(newatomname
, at
->nr
+nadd
);
1909 srenew(newvsite_type
, at
->nr
+nadd
);
1910 srenew(newcgnr
, at
->nr
+nadd
);
1912 for (j
= 0; j
< NMASS
; j
++)
1914 newatomname
[at
->nr
+nadd
-1-j
] = nullptr;
1917 /* calculate starting position for the masses */
1919 /* get atom masses, and set Heavy and Hatoms mass to zero */
1920 for (j
= 0; j
< nrHatoms
; j
++)
1922 mHtot
+= get_amass(Hatoms
[j
], at
, nrtp
, rtp
, rt
);
1923 at
->atom
[Hatoms
[j
]].m
= at
->atom
[Hatoms
[j
]].mB
= 0;
1925 mtot
= mHtot
+ get_amass(Heavy
, at
, nrtp
, rtp
, rt
);
1926 at
->atom
[Heavy
].m
= at
->atom
[Heavy
].mB
= 0;
1932 fact
= std::sqrt(fact2
);
1933 /* generate vectors parallel and perpendicular to rotational axis:
1934 * rpar = Heavy -> Hcom
1935 * rperp = Hcom -> H1 */
1937 for (j
= 0; j
< nrHatoms
; j
++)
1939 rvec_inc(rpar
, (*x
)[Hatoms
[j
]]);
1941 svmul(1.0/nrHatoms
, rpar
, rpar
); /* rpar = ( H1+H2+H3 ) / 3 */
1942 rvec_dec(rpar
, (*x
)[Heavy
]); /* - Heavy */
1943 rvec_sub((*x
)[Hatoms
[0]], (*x
)[Heavy
], rperp
);
1944 rvec_dec(rperp
, rpar
); /* rperp = H1 - Heavy - rpar */
1945 /* calc mass positions */
1946 svmul(fact2
, rpar
, temp
);
1947 for (j
= 0; (j
< NMASS
); j
++) /* xM = xN + fact2 * rpar +/- fact * rperp */
1949 rvec_add((*x
)[Heavy
], temp
, newx
[ni0
+j
]);
1951 svmul(fact
, rperp
, temp
);
1952 rvec_inc(newx
[ni0
], temp
);
1953 rvec_dec(newx
[ni0
+1], temp
);
1954 /* set atom parameters for the masses */
1955 for (j
= 0; (j
< NMASS
); j
++)
1957 /* make name: "M??#" or "M?#" (? is atomname, # is number) */
1959 for (k
= 0; (*at
->atomname
[Heavy
])[k
] && ( k
< NMASS
); k
++)
1961 name
[k
+1] = (*at
->atomname
[Heavy
])[k
];
1963 name
[k
+1] = atomnamesuffix
[j
];
1965 newatomname
[ni0
+j
] = put_symtab(symtab
, name
);
1966 newatom
[ni0
+j
].m
= newatom
[ni0
+j
].mB
= mtot
/NMASS
;
1967 newatom
[ni0
+j
].q
= newatom
[ni0
+j
].qB
= 0.0;
1968 newatom
[ni0
+j
].type
= newatom
[ni0
+j
].typeB
= tpM
;
1969 newatom
[ni0
+j
].ptype
= eptAtom
;
1970 newatom
[ni0
+j
].resind
= at
->atom
[i0
].resind
;
1971 newatom
[ni0
+j
].elem
[0] = 'M';
1972 newatom
[ni0
+j
].elem
[1] = '\0';
1973 newvsite_type
[ni0
+j
] = NOTSET
;
1974 newcgnr
[ni0
+j
] = (*cgnr
)[i0
];
1976 /* add constraints between dummy masses and to heavies[0] */
1977 /* 'add_shift' says which atoms won't be renumbered afterwards */
1978 my_add_param(&(plist
[F_CONSTRNC
]), heavies
[0], add_shift
+ni0
, NOTSET
);
1979 my_add_param(&(plist
[F_CONSTRNC
]), heavies
[0], add_shift
+ni0
+1, NOTSET
);
1980 my_add_param(&(plist
[F_CONSTRNC
]), add_shift
+ni0
, add_shift
+ni0
+1, NOTSET
);
1982 /* generate Heavy, H1, H2 and H3 from M1, M2 and heavies[0] */
1983 /* note that vsite_type cannot be NOTSET, because we just set it */
1984 add_vsite3_atoms (&plist
[(*vsite_type
)[Heavy
]],
1985 Heavy
, heavies
[0], add_shift
+ni0
, add_shift
+ni0
+1,
1987 for (j
= 0; j
< nrHatoms
; j
++)
1989 add_vsite3_atoms(&plist
[(*vsite_type
)[Hatoms
[j
]]],
1990 Hatoms
[j
], heavies
[0], add_shift
+ni0
, add_shift
+ni0
+1,
2004 "Warning: cannot convert atom %d %s (bound to a heavy atom "
2006 " %d bonds and %d bound hydrogens atoms) to virtual site\n",
2007 i
+1, *(at
->atomname
[i
]), tpname
, nrbonds
, nrHatoms
);
2011 /* add vsite parameters to topology,
2012 also get rid of negative vsite_types */
2013 add_vsites(plist
, (*vsite_type
), Heavy
, nrHatoms
, Hatoms
,
2014 nrheavies
, heavies
);
2015 /* transfer mass of virtual site to Heavy atom */
2016 for (j
= 0; j
< nrHatoms
; j
++)
2018 if (is_vsite((*vsite_type
)[Hatoms
[j
]]))
2020 at
->atom
[Heavy
].m
+= at
->atom
[Hatoms
[j
]].m
;
2021 at
->atom
[Heavy
].mB
= at
->atom
[Heavy
].m
;
2022 at
->atom
[Hatoms
[j
]].m
= at
->atom
[Hatoms
[j
]].mB
= 0;
2029 fprintf(debug
, "atom %d: ", o2n
[i
]+1);
2030 print_bonds(debug
, o2n
, nrHatoms
, Hatoms
, Heavy
, nrheavies
, heavies
);
2032 } /* if vsite NOTSET & is hydrogen */
2034 } /* for i < at->nr */
2036 gmx_residuetype_destroy(rt
);
2040 fprintf(debug
, "Before inserting new atoms:\n");
2041 for (i
= 0; i
< at
->nr
; i
++)
2043 fprintf(debug
, "%4d %4d %4s %4d %4s %6d %-10s\n", i
+1, o2n
[i
]+1,
2044 at
->atomname
[i
] ? *(at
->atomname
[i
]) : "(NULL)",
2045 at
->resinfo
[at
->atom
[i
].resind
].nr
,
2046 at
->resinfo
[at
->atom
[i
].resind
].name
?
2047 *(at
->resinfo
[at
->atom
[i
].resind
].name
) : "(NULL)",
2049 ((*vsite_type
)[i
] == NOTSET
) ?
2050 "NOTSET" : interaction_function
[(*vsite_type
)[i
]].name
);
2052 fprintf(debug
, "new atoms to be inserted:\n");
2053 for (i
= 0; i
< at
->nr
+nadd
; i
++)
2057 fprintf(debug
, "%4d %4s %4d %6d %-10s\n", i
+1,
2058 newatomname
[i
] ? *(newatomname
[i
]) : "(NULL)",
2059 newatom
[i
].resind
, newcgnr
[i
],
2060 (newvsite_type
[i
] == NOTSET
) ?
2061 "NOTSET" : interaction_function
[newvsite_type
[i
]].name
);
2066 /* add all original atoms to the new arrays, using o2n index array */
2067 for (i
= 0; i
< at
->nr
; i
++)
2069 newatomname
[o2n
[i
]] = at
->atomname
[i
];
2070 newatom
[o2n
[i
]] = at
->atom
[i
];
2071 newvsite_type
[o2n
[i
]] = (*vsite_type
)[i
];
2072 newcgnr
[o2n
[i
]] = (*cgnr
) [i
];
2073 copy_rvec((*x
)[i
], newx
[o2n
[i
]]);
2075 /* throw away old atoms */
2077 sfree(at
->atomname
);
2081 /* put in the new ones */
2084 at
->atomname
= newatomname
;
2085 *vsite_type
= newvsite_type
;
2088 if (at
->nr
> add_shift
)
2090 gmx_fatal(FARGS
, "Added impossible amount of dummy masses "
2091 "(%d on a total of %d atoms)\n", nadd
, at
->nr
-nadd
);
2096 fprintf(debug
, "After inserting new atoms:\n");
2097 for (i
= 0; i
< at
->nr
; i
++)
2099 fprintf(debug
, "%4d %4s %4d %4s %6d %-10s\n", i
+1,
2100 at
->atomname
[i
] ? *(at
->atomname
[i
]) : "(NULL)",
2101 at
->resinfo
[at
->atom
[i
].resind
].nr
,
2102 at
->resinfo
[at
->atom
[i
].resind
].name
?
2103 *(at
->resinfo
[at
->atom
[i
].resind
].name
) : "(NULL)",
2105 ((*vsite_type
)[i
] == NOTSET
) ?
2106 "NOTSET" : interaction_function
[(*vsite_type
)[i
]].name
);
2110 /* now renumber all the interactions because of the added atoms */
2111 for (ftype
= 0; ftype
< F_NRE
; ftype
++)
2113 params
= &(plist
[ftype
]);
2116 fprintf(debug
, "Renumbering %d %s\n", params
->nr
,
2117 interaction_function
[ftype
].longname
);
2119 for (i
= 0; i
< params
->nr
; i
++)
2121 for (j
= 0; j
< NRAL(ftype
); j
++)
2123 if (params
->param
[i
].a
[j
] >= add_shift
)
2127 fprintf(debug
, " [%d -> %d]", params
->param
[i
].a
[j
],
2128 params
->param
[i
].a
[j
]-add_shift
);
2130 params
->param
[i
].a
[j
] = params
->param
[i
].a
[j
]-add_shift
;
2136 fprintf(debug
, " [%d -> %d]", params
->param
[i
].a
[j
],
2137 o2n
[params
->param
[i
].a
[j
]]);
2139 params
->param
[i
].a
[j
] = o2n
[params
->param
[i
].a
[j
]];
2144 fprintf(debug
, "\n");
2148 /* now check if atoms in the added constraints are in increasing order */
2149 params
= &(plist
[F_CONSTRNC
]);
2150 for (i
= 0; i
< params
->nr
; i
++)
2152 if (params
->param
[i
].ai() > params
->param
[i
].aj())
2154 j
= params
->param
[i
].aj();
2155 params
->param
[i
].aj() = params
->param
[i
].ai();
2156 params
->param
[i
].ai() = j
;
2163 /* tell the user what we did */
2164 fprintf(stderr
, "Marked %d virtual sites\n", nvsite
);
2165 fprintf(stderr
, "Added %d dummy masses\n", nadd
);
2166 fprintf(stderr
, "Added %d new constraints\n", plist
[F_CONSTRNC
].nr
);
2169 void do_h_mass(t_params
*psb
, int vsite_type
[], t_atoms
*at
, real mHmult
,
2170 gmx_bool bDeuterate
)
2174 /* loop over all atoms */
2175 for (i
= 0; i
< at
->nr
; i
++)
2177 /* adjust masses if i is hydrogen and not a virtual site */
2178 if (!is_vsite(vsite_type
[i
]) && is_hydrogen(*(at
->atomname
[i
])) )
2180 /* find bonded heavy atom */
2182 for (j
= 0; (j
< psb
->nr
) && (a
== NOTSET
); j
++)
2184 /* if other atom is not a virtual site, it is the one we want */
2185 if ( (psb
->param
[j
].ai() == i
) &&
2186 !is_vsite(vsite_type
[psb
->param
[j
].aj()]) )
2188 a
= psb
->param
[j
].aj();
2190 else if ( (psb
->param
[j
].aj() == i
) &&
2191 !is_vsite(vsite_type
[psb
->param
[j
].ai()]) )
2193 a
= psb
->param
[j
].ai();
2198 gmx_fatal(FARGS
, "Unbound hydrogen atom (%d) found while adjusting mass",
2202 /* adjust mass of i (hydrogen) with mHmult
2203 and correct mass of a (bonded atom) with same amount */
2206 at
->atom
[a
].m
-= (mHmult
-1.0)*at
->atom
[i
].m
;
2207 at
->atom
[a
].mB
-= (mHmult
-1.0)*at
->atom
[i
].m
;
2209 at
->atom
[i
].m
*= mHmult
;
2210 at
->atom
[i
].mB
*= mHmult
;