?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Replaced NOTSET from typedefs.h by local solutions.
[gromacs.git] / src / gromacs / gmxpreprocess / gen_vsite.cpp
blobd7164b6b715ec90064580fed85c02d2b956dc6a9
1 /*
2 * This file is part of the GROMACS molecular simulation package.
3 *
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, 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.
10 *
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.
15 *
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.
20 *
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.
25 *
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.
33 *
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.
36 */
37 #include "gmxpre.h"
38
39 #include "gen_vsite.h"
40
41 #include <stdio.h>
42 #include <stdlib.h>
43 #include <string.h>
44
45 #include <cmath>
46
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/legacyheaders/names.h"
55 #include "gromacs/legacyheaders/types/ifunc.h"
56 #include "gromacs/math/units.h"
57 #include "gromacs/math/vec.h"
58 #include "gromacs/topology/residuetypes.h"
59 #include "gromacs/topology/symtab.h"
60 #include "gromacs/utility/basedefinitions.h"
61 #include "gromacs/utility/cstringutil.h"
62 #include "gromacs/utility/fatalerror.h"
63 #include "gromacs/utility/futil.h"
64 #include "gromacs/utility/real.h"
65 #include "gromacs/utility/smalloc.h"
66
67 #define MAXNAME 32
68 #define OPENDIR '[' /* starting sign for directive */
69 #define CLOSEDIR ']' /* ending sign for directive */
70
71 typedef struct {
72 char atomtype[MAXNAME]; /* Type for the XH3/XH2 atom */
73 gmx_bool isplanar; /* If true, the atomtype above and the three connected
74 * ones are in a planar geometry. The two next entries
75 * are undefined in that case
76 */
77 int nhydrogens; /* number of connected hydrogens */
78 char nextheavytype[MAXNAME]; /* Type for the heavy atom bonded to XH2/XH3 */
79 char dummymass[MAXNAME]; /* The type of MNH* or MCH3* dummy mass to use */
80 } t_vsiteconf;
81
82
83 /* Structure to represent average bond and angles values in vsite aromatic
84 * residues. Note that these are NOT necessarily the bonds and angles from the
85 * forcefield; many forcefields (like Amber, OPLS) have some inherent strain in
86 * 5-rings (i.e. the sum of angles is !=540, but impropers keep it planar)
87 */
88 typedef struct {
89 char resname[MAXNAME];
90 int nbonds;
91 int nangles;
92 struct vsitetop_bond {
93 char atom1[MAXNAME];
94 char atom2[MAXNAME];
95 float value;
96 } *bond; /* list of bonds */
97 struct vsitetop_angle {
98 char atom1[MAXNAME];
99 char atom2[MAXNAME];
100 char atom3[MAXNAME];
101 float value;
102 } *angle; /* list of angles */
103 } t_vsitetop;
104
105
106 enum {
107 DDB_CH3, DDB_NH3, DDB_NH2, DDB_PHE, DDB_TYR,
108 DDB_TRP, DDB_HISA, DDB_HISB, DDB_HISH, DDB_DIR_NR
109 };
110
111 typedef char t_dirname[STRLEN];
112
113 static const t_dirname ddb_dirnames[DDB_DIR_NR] = {
114 "CH3",
115 "NH3",
116 "NH2",
117 "PHE",
118 "TYR",
119 "TRP",
120 "HISA",
121 "HISB",
122 "HISH"
123 };
124
125 static int ddb_name2dir(char *name)
126 {
127 /* Translate a directive name to the number of the directive.
128 * HID/HIE/HIP names are translated to the ones we use in Gromacs.
129 */
130
131 int i, index;
132
133 index = -1;
134
135 for (i = 0; i < DDB_DIR_NR && index < 0; i++)
136 {
137 if (!gmx_strcasecmp(name, ddb_dirnames[i]))
138 {
139 index = i;
140 }
141 }
142
143 return index;
144 }
145
146
147 static void read_vsite_database(const char *ddbname,
148 t_vsiteconf **pvsiteconflist, int *nvsiteconf,
149 t_vsitetop **pvsitetoplist, int *nvsitetop)
150 {
151 /* This routine is a quick hack to fix the problem with hardcoded atomtypes
152 * and aromatic vsite parameters by reading them from a ff???.vsd file.
153 *
154 * The file can contain sections [ NH3 ], [ CH3 ], [ NH2 ], and ring residue names.
155 * For the NH3 and CH3 section each line has three fields. The first is the atomtype
156 * (nb: not bonded type) of the N/C atom to be replaced, the second field is
157 * the type of the next heavy atom it is bonded to, and the third field the type
158 * of dummy mass that will be used for this group.
159 *
160 * If the NH2 group planar (sp2 N) a different vsite construct is used, so in this
161 * case the second field should just be the word planar.
162 */
163
164 FILE *ddb;
165 char dirstr[STRLEN];
166 char pline[STRLEN];
167 int i, n, k, nvsite, ntop, curdir;
168 t_vsiteconf *vsiteconflist;
169 t_vsitetop *vsitetoplist;
170 char *ch;
171 char s1[MAXNAME], s2[MAXNAME], s3[MAXNAME], s4[MAXNAME];
172
173 ddb = libopen(ddbname);
174
175 nvsite = *nvsiteconf;
176 vsiteconflist = *pvsiteconflist;
177 ntop = *nvsitetop;
178 vsitetoplist = *pvsitetoplist;
179
180 curdir = -1;
181
182 snew(vsiteconflist, 1);
183 snew(vsitetoplist, 1);
184
185 while (fgets2(pline, STRLEN-2, ddb) != NULL)
186 {
187 strip_comment(pline);
188 trim(pline);
189 if (strlen(pline) > 0)
190 {
191 if (pline[0] == OPENDIR)
192 {
193 strncpy(dirstr, pline+1, STRLEN-2);
194 if ((ch = strchr (dirstr, CLOSEDIR)) != NULL)
195 {
196 (*ch) = 0;
197 }
198 trim (dirstr);
199
200 if (!gmx_strcasecmp(dirstr, "HID") ||
201 !gmx_strcasecmp(dirstr, "HISD"))
202 {
203 sprintf(dirstr, "HISA");
204 }
205 else if (!gmx_strcasecmp(dirstr, "HIE") ||
206 !gmx_strcasecmp(dirstr, "HISE"))
207 {
208 sprintf(dirstr, "HISB");
209 }
210 else if (!gmx_strcasecmp(dirstr, "HIP"))
211 {
212 sprintf(dirstr, "HISH");
213 }
214
215 curdir = ddb_name2dir(dirstr);
216 if (curdir < 0)
217 {
218 gmx_fatal(FARGS, "Invalid directive %s in vsite database %s",
219 dirstr, ddbname);
220 }
221 }
222 else
223 {
224 switch (curdir)
225 {
226 case -1:
227 gmx_fatal(FARGS, "First entry in vsite database must be a directive.\n");
228 break;
229 case DDB_CH3:
230 case DDB_NH3:
231 case DDB_NH2:
232 n = sscanf(pline, "%s%s%s", s1, s2, s3);
233 if (n < 3 && !gmx_strcasecmp(s2, "planar"))
234 {
235 srenew(vsiteconflist, nvsite+1);
236 strncpy(vsiteconflist[nvsite].atomtype, s1, MAXNAME-1);
237 vsiteconflist[nvsite].isplanar = TRUE;
238 vsiteconflist[nvsite].nextheavytype[0] = 0;
239 vsiteconflist[nvsite].dummymass[0] = 0;
240 vsiteconflist[nvsite].nhydrogens = 2;
241 nvsite++;
242 }
243 else if (n == 3)
244 {
245 srenew(vsiteconflist, (nvsite+1));
246 strncpy(vsiteconflist[nvsite].atomtype, s1, MAXNAME-1);
247 vsiteconflist[nvsite].isplanar = FALSE;
248 strncpy(vsiteconflist[nvsite].nextheavytype, s2, MAXNAME-1);
249 strncpy(vsiteconflist[nvsite].dummymass, s3, MAXNAME-1);
250 if (curdir == DDB_NH2)
251 {
252 vsiteconflist[nvsite].nhydrogens = 2;
253 }
254 else
255 {
256 vsiteconflist[nvsite].nhydrogens = 3;
257 }
258 nvsite++;
259 }
260 else
261 {
262 gmx_fatal(FARGS, "Not enough directives in vsite database line: %s\n", pline);
263 }
264 break;
265 case DDB_PHE:
266 case DDB_TYR:
267 case DDB_TRP:
268 case DDB_HISA:
269 case DDB_HISB:
270 case DDB_HISH:
271 i = 0;
272 while ((i < ntop) && gmx_strcasecmp(dirstr, vsitetoplist[i].resname))
273 {
274 i++;
275 }
276 /* Allocate a new topology entry if this is a new residue */
277 if (i == ntop)
278 {
279 srenew(vsitetoplist, ntop+1);
280 ntop++; /* i still points to current vsite topology entry */
281 strncpy(vsitetoplist[i].resname, dirstr, MAXNAME-1);
282 vsitetoplist[i].nbonds = vsitetoplist[i].nangles = 0;
283 snew(vsitetoplist[i].bond, 1);
284 snew(vsitetoplist[i].angle, 1);
285 }
286 n = sscanf(pline, "%s%s%s%s", s1, s2, s3, s4);
287 if (n == 3)
288 {
289 /* bond */
290 k = vsitetoplist[i].nbonds++;
291 srenew(vsitetoplist[i].bond, k+1);
292 strncpy(vsitetoplist[i].bond[k].atom1, s1, MAXNAME-1);
293 strncpy(vsitetoplist[i].bond[k].atom2, s2, MAXNAME-1);
294 vsitetoplist[i].bond[k].value = strtod(s3, NULL);
295 }
296 else if (n == 4)
297 {
298 /* angle */
299 k = vsitetoplist[i].nangles++;
300 srenew(vsitetoplist[i].angle, k+1);
301 strncpy(vsitetoplist[i].angle[k].atom1, s1, MAXNAME-1);
302 strncpy(vsitetoplist[i].angle[k].atom2, s2, MAXNAME-1);
303 strncpy(vsitetoplist[i].angle[k].atom3, s3, MAXNAME-1);
304 vsitetoplist[i].angle[k].value = strtod(s4, NULL);
305 }
306 else
307 {
308 gmx_fatal(FARGS, "Need 3 or 4 values to specify bond/angle values in %s: %s\n", ddbname, pline);
309 }
310 break;
311 default:
312 gmx_fatal(FARGS, "Didnt find a case for directive %s in read_vsite_database\n", dirstr);
313 }
314 }
315 }
316 }
317
318 *pvsiteconflist = vsiteconflist;
319 *pvsitetoplist = vsitetoplist;
320 *nvsiteconf = nvsite;
321 *nvsitetop = ntop;
322
323 gmx_ffclose(ddb);
324 }
325
326 static int nitrogen_is_planar(t_vsiteconf vsiteconflist[], int nvsiteconf, char atomtype[])
327 {
328 /* Return 1 if atomtype exists in database list and is planar, 0 if not,
329 * and -1 if not found.
330 */
331 int i, res;
332 gmx_bool found = FALSE;
333 for (i = 0; i < nvsiteconf && !found; i++)
334 {
335 found = (!gmx_strcasecmp(vsiteconflist[i].atomtype, atomtype) && (vsiteconflist[i].nhydrogens == 2));
336 }
337 if (found)
338 {
339 res = (vsiteconflist[i-1].isplanar == TRUE);
340 }
341 else
342 {
343 res = -1;
344 }
345
346 return res;
347 }
348
349 static char *get_dummymass_name(t_vsiteconf vsiteconflist[], int nvsiteconf, char atom[], char nextheavy[])
350 {
351 /* Return the dummy mass name if found, or NULL if not set in ddb database */
352 int i;
353 gmx_bool found = FALSE;
354 for (i = 0; i < nvsiteconf && !found; i++)
355 {
356 found = (!gmx_strcasecmp(vsiteconflist[i].atomtype, atom) &&
357 !gmx_strcasecmp(vsiteconflist[i].nextheavytype, nextheavy));
358 }
359 if (found)
360 {
361 return vsiteconflist[i-1].dummymass;
362 }
363 else
364 {
365 return NULL;
366 }
367 }
368
369
370
371 static real get_ddb_bond(t_vsitetop *vsitetop, int nvsitetop,
372 const char res[],
373 const char atom1[], const char atom2[])
374 {
375 int i, j;
376
377 i = 0;
378 while (i < nvsitetop && gmx_strcasecmp(res, vsitetop[i].resname))
379 {
380 i++;
381 }
382 if (i == nvsitetop)
383 {
384 gmx_fatal(FARGS, "No vsite information for residue %s found in vsite database.\n", res);
385 }
386 j = 0;
387 while (j < vsitetop[i].nbonds &&
388 ( strcmp(atom1, vsitetop[i].bond[j].atom1) || strcmp(atom2, vsitetop[i].bond[j].atom2)) &&
389 ( strcmp(atom2, vsitetop[i].bond[j].atom1) || strcmp(atom1, vsitetop[i].bond[j].atom2)))
390 {
391 j++;
392 }
393 if (j == vsitetop[i].nbonds)
394 {
395 gmx_fatal(FARGS, "Couldnt find bond %s-%s for residue %s in vsite database.\n", atom1, atom2, res);
396 }
397
398 return vsitetop[i].bond[j].value;
399 }
400
401
402 static real get_ddb_angle(t_vsitetop *vsitetop, int nvsitetop,
403 const char res[], const char atom1[],
404 const char atom2[], const char atom3[])
405 {
406 int i, j;
407
408 i = 0;
409 while (i < nvsitetop && gmx_strcasecmp(res, vsitetop[i].resname))
410 {
411 i++;
412 }
413 if (i == nvsitetop)
414 {
415 gmx_fatal(FARGS, "No vsite information for residue %s found in vsite database.\n", res);
416 }
417 j = 0;
418 while (j < vsitetop[i].nangles &&
419 ( strcmp(atom1, vsitetop[i].angle[j].atom1) ||
420 strcmp(atom2, vsitetop[i].angle[j].atom2) ||
421 strcmp(atom3, vsitetop[i].angle[j].atom3)) &&
422 ( strcmp(atom3, vsitetop[i].angle[j].atom1) ||
423 strcmp(atom2, vsitetop[i].angle[j].atom2) ||
424 strcmp(atom1, vsitetop[i].angle[j].atom3)))
425 {
426 j++;
427 }
428 if (j == vsitetop[i].nangles)
429 {
430 gmx_fatal(FARGS, "Couldnt find angle %s-%s-%s for residue %s in vsite database.\n", atom1, atom2, atom3, res);
431 }
432
433 return vsitetop[i].angle[j].value;
434 }
435
436
437 static void count_bonds(int atom, t_params *psb, char ***atomname,
438 int *nrbonds, int *nrHatoms, int Hatoms[], int *Heavy,
439 int *nrheavies, int heavies[])
440 {
441 int i, heavy, other, nrb, nrH, nrhv;
442
443 /* find heavy atom bound to this hydrogen */
444 heavy = NOTSET;
445 for (i = 0; (i < psb->nr) && (heavy == NOTSET); i++)
446 {
447 if (psb->param[i].ai() == atom)
448 {
449 heavy = psb->param[i].aj();
450 }
451 else if (psb->param[i].aj() == atom)
452 {
453 heavy = psb->param[i].ai();
454 }
455 }
456 if (heavy == NOTSET)
457 {
458 gmx_fatal(FARGS, "unbound hydrogen atom %d", atom+1);
459 }
460 /* find all atoms bound to heavy atom */
461 other = NOTSET;
462 nrb = 0;
463 nrH = 0;
464 nrhv = 0;
465 for (i = 0; i < psb->nr; i++)
466 {
467 if (psb->param[i].ai() == heavy)
468 {
469 other = psb->param[i].aj();
470 }
471 else if (psb->param[i].aj() == heavy)
472 {
473 other = psb->param[i].ai();
474 }
475 if (other != NOTSET)
476 {
477 nrb++;
478 if (is_hydrogen(*(atomname[other])))
479 {
480 Hatoms[nrH] = other;
481 nrH++;
482 }
483 else
484 {
485 heavies[nrhv] = other;
486 nrhv++;
487 }
488 other = NOTSET;
489 }
490 }
491 *Heavy = heavy;
492 *nrbonds = nrb;
493 *nrHatoms = nrH;
494 *nrheavies = nrhv;
495 }
496
497 static void print_bonds(FILE *fp, int o2n[],
498 int nrHatoms, int Hatoms[], int Heavy,
499 int nrheavies, int heavies[])
500 {
501 int i;
502
503 fprintf(fp, "Found: %d Hatoms: ", nrHatoms);
504 for (i = 0; i < nrHatoms; i++)
505 {
506 fprintf(fp, " %d", o2n[Hatoms[i]]+1);
507 }
508 fprintf(fp, "; %d Heavy atoms: %d", nrheavies+1, o2n[Heavy]+1);
509 for (i = 0; i < nrheavies; i++)
510 {
511 fprintf(fp, " %d", o2n[heavies[i]]+1);
512 }
513 fprintf(fp, "\n");
514 }
515
516 static int get_atype(int atom, t_atoms *at, int nrtp, t_restp rtp[],
517 gmx_residuetype_t *rt)
518 {
519 int type;
520 gmx_bool bNterm;
521 int j;
522 t_restp *rtpp;
523
524 if (at->atom[atom].m)
525 {
526 type = at->atom[atom].type;
527 }
528 else
529 {
530 /* get type from rtp */
531 rtpp = get_restp(*(at->resinfo[at->atom[atom].resind].name), nrtp, rtp);
532 bNterm = gmx_residuetype_is_protein(rt, *(at->resinfo[at->atom[atom].resind].name)) &&
533 (at->atom[atom].resind == 0);
534 j = search_jtype(rtpp, *(at->atomname[atom]), bNterm);
535 type = rtpp->atom[j].type;
536 }
537 return type;
538 }
539
540 static int vsite_nm2type(const char *name, gpp_atomtype_t atype)
541 {
542 int tp;
543
544 tp = get_atomtype_type(name, atype);
545 if (tp == NOTSET)
546 {
547 gmx_fatal(FARGS, "Dummy mass type (%s) not found in atom type database",
548 name);
549 }
550
551 return tp;
552 }
553
554 static real get_amass(int atom, t_atoms *at, int nrtp, t_restp rtp[],
555 gmx_residuetype_t *rt)
556 {
557 real mass;
558 gmx_bool bNterm;
559 int j;
560 t_restp *rtpp;
561
562 if (at->atom[atom].m)
563 {
564 mass = at->atom[atom].m;
565 }
566 else
567 {
568 /* get mass from rtp */
569 rtpp = get_restp(*(at->resinfo[at->atom[atom].resind].name), nrtp, rtp);
570 bNterm = gmx_residuetype_is_protein(rt, *(at->resinfo[at->atom[atom].resind].name)) &&
571 (at->atom[atom].resind == 0);
572 j = search_jtype(rtpp, *(at->atomname[atom]), bNterm);
573 mass = rtpp->atom[j].m;
574 }
575 return mass;
576 }
577
578 static void my_add_param(t_params *plist, int ai, int aj, real b)
579 {
580 static real c[MAXFORCEPARAM] =
581 { NOTSET, NOTSET, NOTSET, NOTSET, NOTSET, NOTSET };
582
583 c[0] = b;
584 add_param(plist, ai, aj, c, NULL);
585 }
586
587 static void add_vsites(t_params plist[], int vsite_type[],
588 int Heavy, int nrHatoms, int Hatoms[],
589 int nrheavies, int heavies[])
590 {
591 int i, j, ftype, other, moreheavy;
592 gmx_bool bSwapParity;
593
594 for (i = 0; i < nrHatoms; i++)
595 {
596 ftype = vsite_type[Hatoms[i]];
597 /* Errors in setting the vsite_type should really be caugth earlier,
598 * because here it's not possible to print any useful error message.
599 * But it's still better to print a message than to segfault.
600 */
601 if (ftype == NOTSET)
602 {
603 gmx_incons("Undetected error in setting up virtual sites");
604 }
605 bSwapParity = (ftype < 0);
606 vsite_type[Hatoms[i]] = ftype = abs(ftype);
607 if (ftype == F_BONDS)
608 {
609 if ( (nrheavies != 1) && (nrHatoms != 1) )
610 {
611 gmx_fatal(FARGS, "cannot make constraint in add_vsites for %d heavy "
612 "atoms and %d hydrogen atoms", nrheavies, nrHatoms);
613 }
614 my_add_param(&(plist[F_CONSTRNC]), Hatoms[i], heavies[0], NOTSET);
615 }
616 else
617 {
618 switch (ftype)
619 {
620 case F_VSITE3:
621 case F_VSITE3FD:
622 case F_VSITE3OUT:
623 if (nrheavies < 2)
624 {
625 gmx_fatal(FARGS, "Not enough heavy atoms (%d) for %s (min 3)",
626 nrheavies+1,
627 interaction_function[vsite_type[Hatoms[i]]].name);
628 }
629 add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], heavies[1],
630 bSwapParity);
631 break;
632 case F_VSITE3FAD:
633 {
634 if (nrheavies > 1)
635 {
636 moreheavy = heavies[1];
637 }
638 else
639 {
640 /* find more heavy atoms */
641 other = moreheavy = NOTSET;
642 for (j = 0; (j < plist[F_BONDS].nr) && (moreheavy == NOTSET); j++)
643 {
644 if (plist[F_BONDS].param[j].ai() == heavies[0])
645 {
646 other = plist[F_BONDS].param[j].aj();
647 }
648 else if (plist[F_BONDS].param[j].aj() == heavies[0])
649 {
650 other = plist[F_BONDS].param[j].ai();
651 }
652 if ( (other != NOTSET) && (other != Heavy) )
653 {
654 moreheavy = other;
655 }
656 }
657 if (moreheavy == NOTSET)
658 {
659 gmx_fatal(FARGS, "Unbound molecule part %d-%d", Heavy+1, Hatoms[0]+1);
660 }
661 }
662 add_vsite3_atoms(&plist[ftype], Hatoms[i], Heavy, heavies[0], moreheavy,
663 bSwapParity);
664 break;
665 }
666 case F_VSITE4FD:
667 case F_VSITE4FDN:
668 if (nrheavies < 3)
669 {
670 gmx_fatal(FARGS, "Not enough heavy atoms (%d) for %s (min 4)",
671 nrheavies+1,
672 interaction_function[vsite_type[Hatoms[i]]].name);
673 }
674 add_vsite4_atoms(&plist[ftype],
675 Hatoms[0], Heavy, heavies[0], heavies[1], heavies[2]);
676 break;
677
678 default:
679 gmx_fatal(FARGS, "can't use add_vsites for interaction function %s",
680 interaction_function[vsite_type[Hatoms[i]]].name);
681 } /* switch ftype */
682 } /* else */
683 } /* for i */
684 }
685
686 #define ANGLE_6RING (DEG2RAD*120)
687
688 /* cosine rule: a^2 = b^2 + c^2 - 2 b c cos(alpha) */
689 /* get a^2 when a, b and alpha are given: */
690 #define cosrule(b, c, alpha) ( sqr(b) + sqr(c) - 2*b*c*cos(alpha) )
691 /* get cos(alpha) when a, b and c are given: */
692 #define acosrule(a, b, c) ( (sqr(b)+sqr(c)-sqr(a))/(2*b*c) )
693
694 static int gen_vsites_6ring(t_atoms *at, int *vsite_type[], t_params plist[],
695 int nrfound, int *ats, real bond_cc, real bond_ch,
696 real xcom, gmx_bool bDoZ)
697 {
698 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
699 enum {
700 atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
701 atCZ, atHZ, atNR
702 };
703
704 int i, nvsite;
705 real a, b, dCGCE, tmp1, tmp2, mtot, mG, mrest;
706 real xCG;
707 /* CG, CE1 and CE2 stay and each get a part of the total mass,
708 * so the c-o-m stays the same.
709 */
710
711 if (bDoZ)
712 {
713 if (atNR != nrfound)
714 {
715 gmx_incons("Generating vsites on 6-rings");
716 }
717 }
718
719 /* constraints between CG, CE1 and CE2: */
720 dCGCE = std::sqrt( cosrule(bond_cc, bond_cc, ANGLE_6RING) );
721 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE);
722 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE2], dCGCE);
723 my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atCE2], dCGCE);
724
725 /* rest will be vsite3 */
726 mtot = 0;
727 nvsite = 0;
728 for (i = 0; i < (bDoZ ? atNR : atHZ); i++)
729 {
730 mtot += at->atom[ats[i]].m;
731 if (i != atCG && i != atCE1 && i != atCE2 && (bDoZ || (i != atHZ && i != atCZ) ) )
732 {
733 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
734 (*vsite_type)[ats[i]] = F_VSITE3;
735 nvsite++;
736 }
737 }
738 /* Distribute mass so center-of-mass stays the same.
739 * The center-of-mass in the call is defined with x=0 at
740 * the CE1-CE2 bond and y=0 at the line from CG to the middle of CE1-CE2 bond.
741 */
742 xCG = -bond_cc+bond_cc*cos(ANGLE_6RING);
743
744 mG = at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = xcom*mtot/xCG;
745 mrest = mtot-mG;
746 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB =
747 at->atom[ats[atCE2]].m = at->atom[ats[atCE2]].mB = mrest / 2;
748
749 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
750 tmp1 = dCGCE*sin(ANGLE_6RING*0.5);
751 tmp2 = bond_cc*cos(0.5*ANGLE_6RING) + tmp1;
752 tmp1 *= 2;
753 a = b = -bond_ch / tmp1;
754 /* HE1 and HE2: */
755 add_vsite3_param(&plist[F_VSITE3],
756 ats[atHE1], ats[atCE1], ats[atCE2], ats[atCG], a, b);
757 add_vsite3_param(&plist[F_VSITE3],
758 ats[atHE2], ats[atCE2], ats[atCE1], ats[atCG], a, b);
759 /* CD1, CD2 and CZ: */
760 a = b = tmp2 / tmp1;
761 add_vsite3_param(&plist[F_VSITE3],
762 ats[atCD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
763 add_vsite3_param(&plist[F_VSITE3],
764 ats[atCD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
765 if (bDoZ)
766 {
767 add_vsite3_param(&plist[F_VSITE3],
768 ats[atCZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
769 }
770 /* HD1, HD2 and HZ: */
771 a = b = ( bond_ch + tmp2 ) / tmp1;
772 add_vsite3_param(&plist[F_VSITE3],
773 ats[atHD1], ats[atCE2], ats[atCE1], ats[atCG], a, b);
774 add_vsite3_param(&plist[F_VSITE3],
775 ats[atHD2], ats[atCE1], ats[atCE2], ats[atCG], a, b);
776 if (bDoZ)
777 {
778 add_vsite3_param(&plist[F_VSITE3],
779 ats[atHZ], ats[atCG], ats[atCE1], ats[atCE2], a, b);
780 }
781
782 return nvsite;
783 }
784
785 static int gen_vsites_phe(t_atoms *at, int *vsite_type[], t_params plist[],
786 int nrfound, int *ats, t_vsitetop *vsitetop, int nvsitetop)
787 {
788 real bond_cc, bond_ch;
789 real xcom, mtot;
790 int i;
791 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
792 enum {
793 atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
794 atCZ, atHZ, atNR
795 };
796 real x[atNR];
797 /* Aromatic rings have 6-fold symmetry, so we only need one bond length.
798 * (angle is always 120 degrees).
799 */
800 bond_cc = get_ddb_bond(vsitetop, nvsitetop, "PHE", "CD1", "CE1");
801 bond_ch = get_ddb_bond(vsitetop, nvsitetop, "PHE", "CD1", "HD1");
802
803 x[atCG] = -bond_cc+bond_cc*cos(ANGLE_6RING);
804 x[atCD1] = -bond_cc;
805 x[atHD1] = x[atCD1]+bond_ch*cos(ANGLE_6RING);
806 x[atCE1] = 0;
807 x[atHE1] = x[atCE1]-bond_ch*cos(ANGLE_6RING);
808 x[atCD2] = x[atCD1];
809 x[atHD2] = x[atHD1];
810 x[atCE2] = x[atCE1];
811 x[atHE2] = x[atHE1];
812 x[atCZ] = bond_cc*cos(0.5*ANGLE_6RING);
813 x[atHZ] = x[atCZ]+bond_ch;
814
815 xcom = mtot = 0;
816 for (i = 0; i < atNR; i++)
817 {
818 xcom += x[i]*at->atom[ats[i]].m;
819 mtot += at->atom[ats[i]].m;
820 }
821 xcom /= mtot;
822
823 return gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, TRUE);
824 }
825
826 static void calc_vsite3_param(real xd, real yd, real xi, real yi, real xj, real yj,
827 real xk, real yk, real *a, real *b)
828 {
829 /* determine parameters by solving the equation system, since we know the
830 * virtual site coordinates here.
831 */
832 real dx_ij, dx_ik, dy_ij, dy_ik;
833
834 dx_ij = xj-xi;
835 dy_ij = yj-yi;
836 dx_ik = xk-xi;
837 dy_ik = yk-yi;
838
839 *a = ( (xd-xi)*dy_ik - dx_ik*(yd-yi) ) / (dx_ij*dy_ik - dx_ik*dy_ij);
840 *b = ( yd - yi - (*a)*dy_ij ) / dy_ik;
841 }
842
843
844 static int gen_vsites_trp(gpp_atomtype_t atype, rvec *newx[],
845 t_atom *newatom[], char ***newatomname[],
846 int *o2n[], int *newvsite_type[], int *newcgnr[],
847 t_symtab *symtab, int *nadd, rvec x[], int *cgnr[],
848 t_atoms *at, int *vsite_type[], t_params plist[],
849 int nrfound, int *ats, int add_shift,
850 t_vsitetop *vsitetop, int nvsitetop)
851 {
852 #define NMASS 2
853 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
854 enum {
855 atCB, atCG, atCD1, atHD1, atCD2, atNE1, atHE1, atCE2, atCE3, atHE3,
856 atCZ2, atHZ2, atCZ3, atHZ3, atCH2, atHH2, atNR
857 };
858 /* weights for determining the COM's of both rings (M1 and M2): */
859 real mw[NMASS][atNR] = {
860 { 0, 1, 1, 1, 0.5, 1, 1, 0.5, 0, 0,
861 0, 0, 0, 0, 0, 0 },
862 { 0, 0, 0, 0, 0.5, 0, 0, 0.5, 1, 1,
863 1, 1, 1, 1, 1, 1 }
864 };
865
866 real xi[atNR], yi[atNR];
867 real xcom[NMASS], ycom[NMASS], alpha;
868 real b_CD2_CE2, b_NE1_CE2, b_CG_CD2, b_CH2_HH2, b_CE2_CZ2;
869 real b_NE1_HE1, b_CD2_CE3, b_CE3_CZ3, b_CB_CG;
870 real b_CZ2_CH2, b_CZ2_HZ2, b_CD1_HD1, b_CE3_HE3;
871 real b_CG_CD1, b_CZ3_HZ3;
872 real a_NE1_CE2_CD2, a_CE2_CD2_CG, a_CB_CG_CD2, a_CE2_CD2_CE3;
873 real a_CD2_CG_CD1, a_CE2_CZ2_HZ2, a_CZ2_CH2_HH2;
874 real a_CD2_CE2_CZ2, a_CD2_CE3_CZ3, a_CE3_CZ3_HZ3, a_CG_CD1_HD1;
875 real a_CE2_CZ2_CH2, a_HE1_NE1_CE2, a_CD2_CE3_HE3;
876 int atM[NMASS], tpM, i, i0, j, nvsite;
877 real mM[NMASS], dCBM1, dCBM2, dM1M2;
878 real a, b;
879 rvec r_ij, r_ik, t1, t2;
880 char name[10];
881
882 if (atNR != nrfound)
883 {
884 gmx_incons("atom types in gen_vsites_trp");
885 }
886 /* Get geometry from database */
887 b_CD2_CE2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CD2", "CE2");
888 b_NE1_CE2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "NE1", "CE2");
889 b_CG_CD1 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CG", "CD1");
890 b_CG_CD2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CG", "CD2");
891 b_CB_CG = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CB", "CG");
892 b_CE2_CZ2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CE2", "CZ2");
893 b_CD2_CE3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CD2", "CE3");
894 b_CE3_CZ3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CE3", "CZ3");
895 b_CZ2_CH2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CZ2", "CH2");
896
897 b_CD1_HD1 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CD1", "HD1");
898 b_CZ2_HZ2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CZ2", "HZ2");
899 b_NE1_HE1 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "NE1", "HE1");
900 b_CH2_HH2 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CH2", "HH2");
901 b_CE3_HE3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CE3", "HE3");
902 b_CZ3_HZ3 = get_ddb_bond(vsitetop, nvsitetop, "TRP", "CZ3", "HZ3");
903
904 a_NE1_CE2_CD2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "NE1", "CE2", "CD2");
905 a_CE2_CD2_CG = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CD2", "CG");
906 a_CB_CG_CD2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CB", "CG", "CD2");
907 a_CD2_CG_CD1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CG", "CD1");
908 /*a_CB_CG_CD1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CB", "CG", "CD1"); unused */
909
910 a_CE2_CD2_CE3 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CD2", "CE3");
911 a_CD2_CE2_CZ2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CE2", "CZ2");
912 a_CD2_CE3_CZ3 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CE3", "CZ3");
913 a_CE3_CZ3_HZ3 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE3", "CZ3", "HZ3");
914 a_CZ2_CH2_HH2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CZ2", "CH2", "HH2");
915 a_CE2_CZ2_HZ2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CZ2", "HZ2");
916 a_CE2_CZ2_CH2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CE2", "CZ2", "CH2");
917 a_CG_CD1_HD1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CG", "CD1", "HD1");
918 a_HE1_NE1_CE2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "HE1", "NE1", "CE2");
919 a_CD2_CE3_HE3 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TRP", "CD2", "CE3", "HE3");
920
921 /* Calculate local coordinates.
922 * y-axis (x=0) is the bond CD2-CE2.
923 * x-axis (y=0) is perpendicular to the bond CD2-CE2 and
924 * intersects the middle of the bond.
925 */
926 xi[atCD2] = 0;
927 yi[atCD2] = -0.5*b_CD2_CE2;
928
929 xi[atCE2] = 0;
930 yi[atCE2] = 0.5*b_CD2_CE2;
931
932 xi[atNE1] = -b_NE1_CE2*sin(a_NE1_CE2_CD2);
933 yi[atNE1] = yi[atCE2]-b_NE1_CE2*cos(a_NE1_CE2_CD2);
934
935 xi[atCG] = -b_CG_CD2*sin(a_CE2_CD2_CG);
936 yi[atCG] = yi[atCD2]+b_CG_CD2*cos(a_CE2_CD2_CG);
937
938 alpha = a_CE2_CD2_CG + M_PI - a_CB_CG_CD2;
939 xi[atCB] = xi[atCG]-b_CB_CG*sin(alpha);
940 yi[atCB] = yi[atCG]+b_CB_CG*cos(alpha);
941
942 alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - M_PI;
943 xi[atCD1] = xi[atCG]-b_CG_CD1*sin(alpha);
944 yi[atCD1] = yi[atCG]+b_CG_CD1*cos(alpha);
945
946 xi[atCE3] = b_CD2_CE3*sin(a_CE2_CD2_CE3);
947 yi[atCE3] = yi[atCD2]+b_CD2_CE3*cos(a_CE2_CD2_CE3);
948
949 xi[atCZ2] = b_CE2_CZ2*sin(a_CD2_CE2_CZ2);
950 yi[atCZ2] = yi[atCE2]-b_CE2_CZ2*cos(a_CD2_CE2_CZ2);
951
952 alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - M_PI;
953 xi[atCZ3] = xi[atCE3]+b_CE3_CZ3*sin(alpha);
954 yi[atCZ3] = yi[atCE3]+b_CE3_CZ3*cos(alpha);
955
956 alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - M_PI;
957 xi[atCH2] = xi[atCZ2]+b_CZ2_CH2*sin(alpha);
958 yi[atCH2] = yi[atCZ2]-b_CZ2_CH2*cos(alpha);
959
960 /* hydrogens */
961 alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - a_CG_CD1_HD1;
962 xi[atHD1] = xi[atCD1]-b_CD1_HD1*sin(alpha);
963 yi[atHD1] = yi[atCD1]+b_CD1_HD1*cos(alpha);
964
965 alpha = a_NE1_CE2_CD2 + M_PI - a_HE1_NE1_CE2;
966 xi[atHE1] = xi[atNE1]-b_NE1_HE1*sin(alpha);
967 yi[atHE1] = yi[atNE1]-b_NE1_HE1*cos(alpha);
968
969 alpha = a_CE2_CD2_CE3 + M_PI - a_CD2_CE3_HE3;
970 xi[atHE3] = xi[atCE3]+b_CE3_HE3*sin(alpha);
971 yi[atHE3] = yi[atCE3]+b_CE3_HE3*cos(alpha);
972
973 alpha = a_CD2_CE2_CZ2 + M_PI - a_CE2_CZ2_HZ2;
974 xi[atHZ2] = xi[atCZ2]+b_CZ2_HZ2*sin(alpha);
975 yi[atHZ2] = yi[atCZ2]-b_CZ2_HZ2*cos(alpha);
976
977 alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - a_CZ2_CH2_HH2;
978 xi[atHZ3] = xi[atCZ3]+b_CZ3_HZ3*sin(alpha);
979 yi[atHZ3] = yi[atCZ3]+b_CZ3_HZ3*cos(alpha);
980
981 alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - a_CE3_CZ3_HZ3;
982 xi[atHH2] = xi[atCH2]+b_CH2_HH2*sin(alpha);
983 yi[atHH2] = yi[atCH2]-b_CH2_HH2*cos(alpha);
984
985 /* Calculate masses for each ring and put it on the dummy masses */
986 for (j = 0; j < NMASS; j++)
987 {
988 mM[j] = xcom[j] = ycom[j] = 0;
989 }
990 for (i = 0; i < atNR; i++)
991 {
992 if (i != atCB)
993 {
994 for (j = 0; j < NMASS; j++)
995 {
996 mM[j] += mw[j][i] * at->atom[ats[i]].m;
997 xcom[j] += xi[i] * mw[j][i] * at->atom[ats[i]].m;
998 ycom[j] += yi[i] * mw[j][i] * at->atom[ats[i]].m;
999 }
1000 }
1001 }
1002 for (j = 0; j < NMASS; j++)
1003 {
1004 xcom[j] /= mM[j];
1005 ycom[j] /= mM[j];
1006 }
1007
1008 /* get dummy mass type */
1009 tpM = vsite_nm2type("MW", atype);
1010 /* make space for 2 masses: shift all atoms starting with CB */
1011 i0 = ats[atCB];
1012 for (j = 0; j < NMASS; j++)
1013 {
1014 atM[j] = i0+*nadd+j;
1015 }
1016 if (debug)
1017 {
1018 fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, (*o2n)[i0]+1);
1019 }
1020 *nadd += NMASS;
1021 for (j = i0; j < at->nr; j++)
1022 {
1023 (*o2n)[j] = j+*nadd;
1024 }
1025 srenew(*newx, at->nr+*nadd);
1026 srenew(*newatom, at->nr+*nadd);
1027 srenew(*newatomname, at->nr+*nadd);
1028 srenew(*newvsite_type, at->nr+*nadd);
1029 srenew(*newcgnr, at->nr+*nadd);
1030 for (j = 0; j < NMASS; j++)
1031 {
1032 (*newatomname)[at->nr+*nadd-1-j] = NULL;
1033 }
1034
1035 /* Dummy masses will be placed at the center-of-mass in each ring. */
1036
1037 /* calc initial position for dummy masses in real (non-local) coordinates.
1038 * Cheat by using the routine to calculate virtual site parameters. It is
1039 * much easier when we have the coordinates expressed in terms of
1040 * CB, CG, CD2.
1041 */
1042 rvec_sub(x[ats[atCB]], x[ats[atCG]], r_ij);
1043 rvec_sub(x[ats[atCD2]], x[ats[atCG]], r_ik);
1044 calc_vsite3_param(xcom[0], ycom[0], xi[atCG], yi[atCG], xi[atCB], yi[atCB],
1045 xi[atCD2], yi[atCD2], &a, &b);
1046 svmul(a, r_ij, t1);
1047 svmul(b, r_ik, t2);
1048 rvec_add(t1, t2, t1);
1049 rvec_add(t1, x[ats[atCG]], (*newx)[atM[0]]);
1050
1051 calc_vsite3_param(xcom[1], ycom[1], xi[atCG], yi[atCG], xi[atCB], yi[atCB],
1052 xi[atCD2], yi[atCD2], &a, &b);
1053 svmul(a, r_ij, t1);
1054 svmul(b, r_ik, t2);
1055 rvec_add(t1, t2, t1);
1056 rvec_add(t1, x[ats[atCG]], (*newx)[atM[1]]);
1057
1058 /* set parameters for the masses */
1059 for (j = 0; j < NMASS; j++)
1060 {
1061 sprintf(name, "MW%d", j+1);
1062 (*newatomname) [atM[j]] = put_symtab(symtab, name);
1063 (*newatom) [atM[j]].m = (*newatom)[atM[j]].mB = mM[j];
1064 (*newatom) [atM[j]].q = (*newatom)[atM[j]].qB = 0.0;
1065 (*newatom) [atM[j]].type = (*newatom)[atM[j]].typeB = tpM;
1066 (*newatom) [atM[j]].ptype = eptAtom;
1067 (*newatom) [atM[j]].resind = at->atom[i0].resind;
1068 (*newatom) [atM[j]].elem[0] = 'M';
1069 (*newatom) [atM[j]].elem[1] = '\0';
1070 (*newvsite_type)[atM[j]] = NOTSET;
1071 (*newcgnr) [atM[j]] = (*cgnr)[i0];
1072 }
1073 /* renumber cgnr: */
1074 for (i = i0; i < at->nr; i++)
1075 {
1076 (*cgnr)[i]++;
1077 }
1078
1079 /* constraints between CB, M1 and M2 */
1080 /* 'add_shift' says which atoms won't be renumbered afterwards */
1081 dCBM1 = std::sqrt( sqr(xcom[0]-xi[atCB]) + sqr(ycom[0]-yi[atCB]) );
1082 dM1M2 = std::sqrt( sqr(xcom[0]-xcom[1]) + sqr(ycom[0]-ycom[1]) );
1083 dCBM2 = std::sqrt( sqr(xcom[1]-xi[atCB]) + sqr(ycom[1]-yi[atCB]) );
1084 my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift+atM[0], dCBM1);
1085 my_add_param(&(plist[F_CONSTRNC]), ats[atCB], add_shift+atM[1], dCBM2);
1086 my_add_param(&(plist[F_CONSTRNC]), add_shift+atM[0], add_shift+atM[1], dM1M2);
1087
1088 /* rest will be vsite3 */
1089 nvsite = 0;
1090 for (i = 0; i < atNR; i++)
1091 {
1092 if (i != atCB)
1093 {
1094 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
1095 (*vsite_type)[ats[i]] = F_VSITE3;
1096 nvsite++;
1097 }
1098 }
1099
1100 /* now define all vsites from M1, M2, CB, ie:
1101 r_d = r_M1 + a r_M1_M2 + b r_M1_CB */
1102 for (i = 0; i < atNR; i++)
1103 {
1104 if ( (*vsite_type)[ats[i]] == F_VSITE3)
1105 {
1106 calc_vsite3_param(xi[i], yi[i], xcom[0], ycom[0], xcom[1], ycom[1], xi[atCB], yi[atCB], &a, &b);
1107 add_vsite3_param(&plist[F_VSITE3],
1108 ats[i], add_shift+atM[0], add_shift+atM[1], ats[atCB], a, b);
1109 }
1110 }
1111 return nvsite;
1112 #undef NMASS
1113 }
1114
1115
1116 static int gen_vsites_tyr(gpp_atomtype_t atype, rvec *newx[],
1117 t_atom *newatom[], char ***newatomname[],
1118 int *o2n[], int *newvsite_type[], int *newcgnr[],
1119 t_symtab *symtab, int *nadd, rvec x[], int *cgnr[],
1120 t_atoms *at, int *vsite_type[], t_params plist[],
1121 int nrfound, int *ats, int add_shift,
1122 t_vsitetop *vsitetop, int nvsitetop)
1123 {
1124 int nvsite, i, i0, j, atM, tpM;
1125 real dCGCE, dCEOH, dCGM, tmp1, a, b;
1126 real bond_cc, bond_ch, bond_co, bond_oh, angle_coh;
1127 real xcom, mtot;
1128 real vdist, mM;
1129 rvec r1;
1130 char name[10];
1131
1132 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1133 enum {
1134 atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
1135 atCZ, atOH, atHH, atNR
1136 };
1137 real xi[atNR];
1138 /* CG, CE1, CE2 (as in general 6-ring) and OH and HH stay,
1139 rest gets virtualized.
1140 Now we have two linked triangles with one improper keeping them flat */
1141 if (atNR != nrfound)
1142 {
1143 gmx_incons("Number of atom types in gen_vsites_tyr");
1144 }
1145
1146 /* Aromatic rings have 6-fold symmetry, so we only need one bond length
1147 * for the ring part (angle is always 120 degrees).
1148 */
1149 bond_cc = get_ddb_bond(vsitetop, nvsitetop, "TYR", "CD1", "CE1");
1150 bond_ch = get_ddb_bond(vsitetop, nvsitetop, "TYR", "CD1", "HD1");
1151 bond_co = get_ddb_bond(vsitetop, nvsitetop, "TYR", "CZ", "OH");
1152 bond_oh = get_ddb_bond(vsitetop, nvsitetop, "TYR", "OH", "HH");
1153 angle_coh = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, "TYR", "CZ", "OH", "HH");
1154
1155 xi[atCG] = -bond_cc+bond_cc*cos(ANGLE_6RING);
1156 xi[atCD1] = -bond_cc;
1157 xi[atHD1] = xi[atCD1]+bond_ch*cos(ANGLE_6RING);
1158 xi[atCE1] = 0;
1159 xi[atHE1] = xi[atCE1]-bond_ch*cos(ANGLE_6RING);
1160 xi[atCD2] = xi[atCD1];
1161 xi[atHD2] = xi[atHD1];
1162 xi[atCE2] = xi[atCE1];
1163 xi[atHE2] = xi[atHE1];
1164 xi[atCZ] = bond_cc*cos(0.5*ANGLE_6RING);
1165 xi[atOH] = xi[atCZ]+bond_co;
1166
1167 xcom = mtot = 0;
1168 for (i = 0; i < atOH; i++)
1169 {
1170 xcom += xi[i]*at->atom[ats[i]].m;
1171 mtot += at->atom[ats[i]].m;
1172 }
1173 xcom /= mtot;
1174
1175 /* first do 6 ring as default,
1176 except CZ (we'll do that different) and HZ (we don't have that): */
1177 nvsite = gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, FALSE);
1178
1179 /* then construct CZ from the 2nd triangle */
1180 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
1181 a = b = 0.5 * bond_co / ( bond_co - bond_cc*cos(ANGLE_6RING) );
1182 add_vsite3_param(&plist[F_VSITE3],
1183 ats[atCZ], ats[atOH], ats[atCE1], ats[atCE2], a, b);
1184 at->atom[ats[atCZ]].m = at->atom[ats[atCZ]].mB = 0;
1185
1186 /* constraints between CE1, CE2 and OH */
1187 dCGCE = std::sqrt( cosrule(bond_cc, bond_cc, ANGLE_6RING) );
1188 dCEOH = std::sqrt( cosrule(bond_cc, bond_co, ANGLE_6RING) );
1189 my_add_param(&(plist[F_CONSTRNC]), ats[atCE1], ats[atOH], dCEOH);
1190 my_add_param(&(plist[F_CONSTRNC]), ats[atCE2], ats[atOH], dCEOH);
1191
1192 /* We also want to constrain the angle C-O-H, but since CZ is constructed
1193 * we need to introduce a constraint to CG.
1194 * CG is much further away, so that will lead to instabilities in LINCS
1195 * when we constrain both CG-HH and OH-HH distances. Instead of requiring
1196 * the use of lincs_order=8 we introduce a dummy mass three times further
1197 * away from OH than HH. The mass is accordingly a third, with the remaining
1198 * 2/3 moved to OH. This shouldnt cause any problems since the forces will
1199 * apply to the HH constructed atom and not directly on the virtual mass.
1200 */
1201
1202 vdist = 2.0*bond_oh;
1203 mM = at->atom[ats[atHH]].m/2.0;
1204 at->atom[ats[atOH]].m += mM; /* add 1/2 of original H mass */
1205 at->atom[ats[atOH]].mB += mM; /* add 1/2 of original H mass */
1206 at->atom[ats[atHH]].m = at->atom[ats[atHH]].mB = 0;
1207
1208 /* get dummy mass type */
1209 tpM = vsite_nm2type("MW", atype);
1210 /* make space for 1 mass: shift HH only */
1211 i0 = ats[atHH];
1212 atM = i0+*nadd;
1213 if (debug)
1214 {
1215 fprintf(stderr, "Inserting 1 dummy mass at %d\n", (*o2n)[i0]+1);
1216 }
1217 (*nadd)++;
1218 for (j = i0; j < at->nr; j++)
1219 {
1220 (*o2n)[j] = j+*nadd;
1221 }
1222 srenew(*newx, at->nr+*nadd);
1223 srenew(*newatom, at->nr+*nadd);
1224 srenew(*newatomname, at->nr+*nadd);
1225 srenew(*newvsite_type, at->nr+*nadd);
1226 srenew(*newcgnr, at->nr+*nadd);
1227 (*newatomname)[at->nr+*nadd-1] = NULL;
1228
1229 /* Calc the dummy mass initial position */
1230 rvec_sub(x[ats[atHH]], x[ats[atOH]], r1);
1231 svmul(2.0, r1, r1);
1232 rvec_add(r1, x[ats[atHH]], (*newx)[atM]);
1233
1234 strcpy(name, "MW1");
1235 (*newatomname) [atM] = put_symtab(symtab, name);
1236 (*newatom) [atM].m = (*newatom)[atM].mB = mM;
1237 (*newatom) [atM].q = (*newatom)[atM].qB = 0.0;
1238 (*newatom) [atM].type = (*newatom)[atM].typeB = tpM;
1239 (*newatom) [atM].ptype = eptAtom;
1240 (*newatom) [atM].resind = at->atom[i0].resind;
1241 (*newatom) [atM].elem[0] = 'M';
1242 (*newatom) [atM].elem[1] = '\0';
1243 (*newvsite_type)[atM] = NOTSET;
1244 (*newcgnr) [atM] = (*cgnr)[i0];
1245 /* renumber cgnr: */
1246 for (i = i0; i < at->nr; i++)
1247 {
1248 (*cgnr)[i]++;
1249 }
1250
1251 (*vsite_type)[ats[atHH]] = F_VSITE2;
1252 nvsite++;
1253 /* assume we also want the COH angle constrained: */
1254 tmp1 = bond_cc*cos(0.5*ANGLE_6RING) + dCGCE*sin(ANGLE_6RING*0.5) + bond_co;
1255 dCGM = std::sqrt( cosrule(tmp1, vdist, angle_coh) );
1256 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], add_shift+atM, dCGM);
1257 my_add_param(&(plist[F_CONSTRNC]), ats[atOH], add_shift+atM, vdist);
1258
1259 add_vsite2_param(&plist[F_VSITE2],
1260 ats[atHH], ats[atOH], add_shift+atM, 1.0/2.0);
1261 return nvsite;
1262 }
1263
1264 static int gen_vsites_his(t_atoms *at, int *vsite_type[], t_params plist[],
1265 int nrfound, int *ats, t_vsitetop *vsitetop, int nvsitetop)
1266 {
1267 int nvsite, i;
1268 real a, b, alpha, dCGCE1, dCGNE2;
1269 real sinalpha, cosalpha;
1270 real xcom, ycom, mtot;
1271 real mG, mrest, mCE1, mNE2;
1272 real b_CG_ND1, b_ND1_CE1, b_CE1_NE2, b_CG_CD2, b_CD2_NE2;
1273 real b_ND1_HD1, b_NE2_HE2, b_CE1_HE1, b_CD2_HD2;
1274 real a_CG_ND1_CE1, a_CG_CD2_NE2, a_ND1_CE1_NE2, a_CE1_NE2_CD2;
1275 real a_NE2_CE1_HE1, a_NE2_CD2_HD2, a_CE1_ND1_HD1, a_CE1_NE2_HE2;
1276 char resname[10];
1277
1278 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1279 enum {
1280 atCG, atND1, atHD1, atCD2, atHD2, atCE1, atHE1, atNE2, atHE2, atNR
1281 };
1282 real x[atNR], y[atNR];
1283
1284 /* CG, CE1 and NE2 stay, each gets part of the total mass,
1285 rest gets virtualized */
1286 /* check number of atoms, 3 hydrogens may be missing: */
1287 /* assert( nrfound >= atNR-3 || nrfound <= atNR );
1288 * Don't understand the above logic. Shouldn't it be && rather than || ???
1289 */
1290 if ((nrfound < atNR-3) || (nrfound > atNR))
1291 {
1292 gmx_incons("Generating vsites for HIS");
1293 }
1294
1295 /* avoid warnings about uninitialized variables */
1296 b_ND1_HD1 = b_NE2_HE2 = b_CE1_HE1 = b_CD2_HD2 = a_NE2_CE1_HE1 =
1297 a_NE2_CD2_HD2 = a_CE1_ND1_HD1 = a_CE1_NE2_HE2 = 0;
1298
1299 if (ats[atHD1] != NOTSET)
1300 {
1301 if (ats[atHE2] != NOTSET)
1302 {
1303 sprintf(resname, "HISH");
1304 }
1305 else
1306 {
1307 sprintf(resname, "HISA");
1308 }
1309 }
1310 else
1311 {
1312 sprintf(resname, "HISB");
1313 }
1314
1315 /* Get geometry from database */
1316 b_CG_ND1 = get_ddb_bond(vsitetop, nvsitetop, resname, "CG", "ND1");
1317 b_ND1_CE1 = get_ddb_bond(vsitetop, nvsitetop, resname, "ND1", "CE1");
1318 b_CE1_NE2 = get_ddb_bond(vsitetop, nvsitetop, resname, "CE1", "NE2");
1319 b_CG_CD2 = get_ddb_bond(vsitetop, nvsitetop, resname, "CG", "CD2");
1320 b_CD2_NE2 = get_ddb_bond(vsitetop, nvsitetop, resname, "CD2", "NE2");
1321 a_CG_ND1_CE1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CG", "ND1", "CE1");
1322 a_CG_CD2_NE2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CG", "CD2", "NE2");
1323 a_ND1_CE1_NE2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "ND1", "CE1", "NE2");
1324 a_CE1_NE2_CD2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CE1", "NE2", "CD2");
1325
1326 if (ats[atHD1] != NOTSET)
1327 {
1328 b_ND1_HD1 = get_ddb_bond(vsitetop, nvsitetop, resname, "ND1", "HD1");
1329 a_CE1_ND1_HD1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CE1", "ND1", "HD1");
1330 }
1331 if (ats[atHE2] != NOTSET)
1332 {
1333 b_NE2_HE2 = get_ddb_bond(vsitetop, nvsitetop, resname, "NE2", "HE2");
1334 a_CE1_NE2_HE2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "CE1", "NE2", "HE2");
1335 }
1336 if (ats[atHD2] != NOTSET)
1337 {
1338 b_CD2_HD2 = get_ddb_bond(vsitetop, nvsitetop, resname, "CD2", "HD2");
1339 a_NE2_CD2_HD2 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "NE2", "CD2", "HD2");
1340 }
1341 if (ats[atHE1] != NOTSET)
1342 {
1343 b_CE1_HE1 = get_ddb_bond(vsitetop, nvsitetop, resname, "CE1", "HE1");
1344 a_NE2_CE1_HE1 = DEG2RAD*get_ddb_angle(vsitetop, nvsitetop, resname, "NE2", "CE1", "HE1");
1345 }
1346
1347 /* constraints between CG, CE1 and NE1 */
1348 dCGCE1 = std::sqrt( cosrule(b_CG_ND1, b_ND1_CE1, a_CG_ND1_CE1) );
1349 dCGNE2 = std::sqrt( cosrule(b_CG_CD2, b_CD2_NE2, a_CG_CD2_NE2) );
1350
1351 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atCE1], dCGCE1);
1352 my_add_param(&(plist[F_CONSTRNC]), ats[atCG], ats[atNE2], dCGNE2);
1353 /* we already have a constraint CE1-NE2, so we don't add it again */
1354
1355 /* calculate the positions in a local frame of reference.
1356 * The x-axis is the line from CG that makes a right angle
1357 * with the bond CE1-NE2, and the y-axis the bond CE1-NE2.
1358 */
1359 /* First calculate the x-axis intersection with y-axis (=yCE1).
1360 * Get cos(angle CG-CE1-NE2) :
1361 */
1362 cosalpha = acosrule(dCGNE2, dCGCE1, b_CE1_NE2);
1363 x[atCE1] = 0;
1364 y[atCE1] = cosalpha*dCGCE1;
1365 x[atNE2] = 0;
1366 y[atNE2] = y[atCE1]-b_CE1_NE2;
1367 sinalpha = std::sqrt(1-cosalpha*cosalpha);
1368 x[atCG] = -sinalpha*dCGCE1;
1369 y[atCG] = 0;
1370 x[atHE1] = x[atHE2] = x[atHD1] = x[atHD2] = 0;
1371 y[atHE1] = y[atHE2] = y[atHD1] = y[atHD2] = 0;
1372
1373 /* calculate ND1 and CD2 positions from CE1 and NE2 */
1374
1375 x[atND1] = -b_ND1_CE1*sin(a_ND1_CE1_NE2);
1376 y[atND1] = y[atCE1]-b_ND1_CE1*cos(a_ND1_CE1_NE2);
1377
1378 x[atCD2] = -b_CD2_NE2*sin(a_CE1_NE2_CD2);
1379 y[atCD2] = y[atNE2]+b_CD2_NE2*cos(a_CE1_NE2_CD2);
1380
1381 /* And finally the hydrogen positions */
1382 if (ats[atHE1] != NOTSET)
1383 {
1384 x[atHE1] = x[atCE1] + b_CE1_HE1*sin(a_NE2_CE1_HE1);
1385 y[atHE1] = y[atCE1] - b_CE1_HE1*cos(a_NE2_CE1_HE1);
1386 }
1387 /* HD2 - first get (ccw) angle from (positive) y-axis */
1388 if (ats[atHD2] != NOTSET)
1389 {
1390 alpha = a_CE1_NE2_CD2 + M_PI - a_NE2_CD2_HD2;
1391 x[atHD2] = x[atCD2] - b_CD2_HD2*sin(alpha);
1392 y[atHD2] = y[atCD2] + b_CD2_HD2*cos(alpha);
1393 }
1394 if (ats[atHD1] != NOTSET)
1395 {
1396 /* HD1 - first get (cw) angle from (positive) y-axis */
1397 alpha = a_ND1_CE1_NE2 + M_PI - a_CE1_ND1_HD1;
1398 x[atHD1] = x[atND1] - b_ND1_HD1*sin(alpha);
1399 y[atHD1] = y[atND1] - b_ND1_HD1*cos(alpha);
1400 }
1401 if (ats[atHE2] != NOTSET)
1402 {
1403 x[atHE2] = x[atNE2] + b_NE2_HE2*sin(a_CE1_NE2_HE2);
1404 y[atHE2] = y[atNE2] + b_NE2_HE2*cos(a_CE1_NE2_HE2);
1405 }
1406 /* Have all coordinates now */
1407
1408 /* calc center-of-mass; keep atoms CG, CE1, NE2 and
1409 * set the rest to vsite3
1410 */
1411 mtot = xcom = ycom = 0;
1412 nvsite = 0;
1413 for (i = 0; i < atNR; i++)
1414 {
1415 if (ats[i] != NOTSET)
1416 {
1417 mtot += at->atom[ats[i]].m;
1418 xcom += x[i]*at->atom[ats[i]].m;
1419 ycom += y[i]*at->atom[ats[i]].m;
1420 if (i != atCG && i != atCE1 && i != atNE2)
1421 {
1422 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
1423 (*vsite_type)[ats[i]] = F_VSITE3;
1424 nvsite++;
1425 }
1426 }
1427 }
1428 if (nvsite+3 != nrfound)
1429 {
1430 gmx_incons("Generating vsites for HIS");
1431 }
1432
1433 xcom /= mtot;
1434 ycom /= mtot;
1435
1436 /* distribute mass so that com stays the same */
1437 mG = xcom*mtot/x[atCG];
1438 mrest = mtot-mG;
1439 mCE1 = (ycom-y[atNE2])*mrest/(y[atCE1]-y[atNE2]);
1440 mNE2 = mrest-mCE1;
1441
1442 at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = mG;
1443 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB = mCE1;
1444 at->atom[ats[atNE2]].m = at->atom[ats[atNE2]].mB = mNE2;
1445
1446 /* HE1 */
1447 if (ats[atHE1] != NOTSET)
1448 {
1449 calc_vsite3_param(x[atHE1], y[atHE1], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1450 x[atCG], y[atCG], &a, &b);
1451 add_vsite3_param(&plist[F_VSITE3],
1452 ats[atHE1], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1453 }
1454 /* HE2 */
1455 if (ats[atHE2] != NOTSET)
1456 {
1457 calc_vsite3_param(x[atHE2], y[atHE2], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1458 x[atCG], y[atCG], &a, &b);
1459 add_vsite3_param(&plist[F_VSITE3],
1460 ats[atHE2], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1461 }
1462
1463 /* ND1 */
1464 calc_vsite3_param(x[atND1], y[atND1], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1465 x[atCG], y[atCG], &a, &b);
1466 add_vsite3_param(&plist[F_VSITE3],
1467 ats[atND1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1468
1469 /* CD2 */
1470 calc_vsite3_param(x[atCD2], y[atCD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1471 x[atCG], y[atCG], &a, &b);
1472 add_vsite3_param(&plist[F_VSITE3],
1473 ats[atCD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1474
1475 /* HD1 */
1476 if (ats[atHD1] != NOTSET)
1477 {
1478 calc_vsite3_param(x[atHD1], y[atHD1], x[atNE2], y[atNE2], x[atCE1], y[atCE1],
1479 x[atCG], y[atCG], &a, &b);
1480 add_vsite3_param(&plist[F_VSITE3],
1481 ats[atHD1], ats[atNE2], ats[atCE1], ats[atCG], a, b);
1482 }
1483 /* HD2 */
1484 if (ats[atHD2] != NOTSET)
1485 {
1486 calc_vsite3_param(x[atHD2], y[atHD2], x[atCE1], y[atCE1], x[atNE2], y[atNE2],
1487 x[atCG], y[atCG], &a, &b);
1488 add_vsite3_param(&plist[F_VSITE3],
1489 ats[atHD2], ats[atCE1], ats[atNE2], ats[atCG], a, b);
1490 }
1491 return nvsite;
1492 }
1493
1494 static gmx_bool is_vsite(int vsite_type)
1495 {
1496 if (vsite_type == NOTSET)
1497 {
1498 return FALSE;
1499 }
1500 switch (abs(vsite_type) )
1501 {
1502 case F_VSITE3:
1503 case F_VSITE3FD:
1504 case F_VSITE3OUT:
1505 case F_VSITE3FAD:
1506 case F_VSITE4FD:
1507 case F_VSITE4FDN:
1508 return TRUE;
1509 default:
1510 return FALSE;
1511 }
1512 }
1513
1514 static char atomnamesuffix[] = "1234";
1515
1516 void do_vsites(int nrtp, t_restp rtp[], gpp_atomtype_t atype,
1517 t_atoms *at, t_symtab *symtab, rvec *x[],
1518 t_params plist[], int *vsite_type[], int *cgnr[],
1519 real mHmult, gmx_bool bVsiteAromatics,
1520 const char *ffdir)
1521 {
1522 #define MAXATOMSPERRESIDUE 16
1523 int i, j, k, m, i0, ni0, whatres, resind, add_shift, ftype, nvsite, nadd;
1524 int ai, aj, ak, al;
1525 int nrfound = 0, needed, nrbonds, nrHatoms, Heavy, nrheavies, tpM, tpHeavy;
1526 int Hatoms[4], heavies[4];
1527 gmx_bool bWARNING, bAddVsiteParam, bFirstWater;
1528 matrix tmpmat;
1529 gmx_bool *bResProcessed;
1530 real mHtot, mtot, fact, fact2;
1531 rvec rpar, rperp, temp;
1532 char name[10], tpname[32], nexttpname[32], *ch;
1533 rvec *newx;
1534 int *o2n, *newvsite_type, *newcgnr, ats[MAXATOMSPERRESIDUE];
1535 t_atom *newatom;
1536 t_params *params;
1537 char ***newatomname;
1538 char *resnm = NULL;
1539 int ndb, f;
1540 char **db;
1541 int nvsiteconf, nvsitetop, cmplength;
1542 gmx_bool isN, planarN, bFound;
1543 gmx_residuetype_t*rt;
1544
1545 t_vsiteconf *vsiteconflist;
1546 /* pointer to a list of CH3/NH3/NH2 configuration entries.
1547 * See comments in read_vsite_database. It isnt beautiful,
1548 * but it had to be fixed, and I dont even want to try to
1549 * maintain this part of the code...
1550 */
1551 t_vsitetop *vsitetop;
1552 /* Pointer to a list of geometry (bond/angle) entries for
1553 * residues like PHE, TRP, TYR, HIS, etc., where we need
1554 * to know the geometry to construct vsite aromatics.
1555 * Note that equilibrium geometry isnt necessarily the same
1556 * as the individual bond and angle values given in the
1557 * force field (rings can be strained).
1558 */
1559
1560 /* if bVsiteAromatics=TRUE do_vsites will specifically convert atoms in
1561 PHE, TRP, TYR and HIS to a construction of virtual sites */
1562 enum {
1563 resPHE, resTRP, resTYR, resHIS, resNR
1564 };
1565 const char *resnms[resNR] = { "PHE", "TRP", "TYR", "HIS" };
1566 /* Amber03 alternative names for termini */
1567 const char *resnmsN[resNR] = { "NPHE", "NTRP", "NTYR", "NHIS" };
1568 const char *resnmsC[resNR] = { "CPHE", "CTRP", "CTYR", "CHIS" };
1569 /* HIS can be known as HISH, HIS1, HISA, HID, HIE, HIP, etc. too */
1570 gmx_bool bPartial[resNR] = { FALSE, FALSE, FALSE, TRUE };
1571 /* the atnms for every residue MUST correspond to the enums in the
1572 gen_vsites_* (one for each residue) routines! */
1573 /* also the atom names in atnms MUST be in the same order as in the .rtp! */
1574 const char *atnms[resNR][MAXATOMSPERRESIDUE+1] = {
1575 { "CG", /* PHE */
1576 "CD1", "HD1", "CD2", "HD2",
1577 "CE1", "HE1", "CE2", "HE2",
1578 "CZ", "HZ", NULL },
1579 { "CB", /* TRP */
1580 "CG",
1581 "CD1", "HD1", "CD2",
1582 "NE1", "HE1", "CE2", "CE3", "HE3",
1583 "CZ2", "HZ2", "CZ3", "HZ3",
1584 "CH2", "HH2", NULL },
1585 { "CG", /* TYR */
1586 "CD1", "HD1", "CD2", "HD2",
1587 "CE1", "HE1", "CE2", "HE2",
1588 "CZ", "OH", "HH", NULL },
1589 { "CG", /* HIS */
1590 "ND1", "HD1", "CD2", "HD2",
1591 "CE1", "HE1", "NE2", "HE2", NULL }
1592 };
1593
1594 if (debug)
1595 {
1596 printf("Searching for atoms to make virtual sites ...\n");
1597 fprintf(debug, "# # # VSITES # # #\n");
1598 }
1599
1600 ndb = fflib_search_file_end(ffdir, ".vsd", FALSE, &db);
1601 nvsiteconf = 0;
1602 vsiteconflist = NULL;
1603 nvsitetop = 0;
1604 vsitetop = NULL;
1605 for (f = 0; f < ndb; f++)
1606 {
1607 read_vsite_database(db[f], &vsiteconflist, &nvsiteconf, &vsitetop, &nvsitetop);
1608 sfree(db[f]);
1609 }
1610 sfree(db);
1611
1612 bFirstWater = TRUE;
1613 nvsite = 0;
1614 nadd = 0;
1615 /* we need a marker for which atoms should *not* be renumbered afterwards */
1616 add_shift = 10*at->nr;
1617 /* make arrays where masses can be inserted into */
1618 snew(newx, at->nr);
1619 snew(newatom, at->nr);
1620 snew(newatomname, at->nr);
1621 snew(newvsite_type, at->nr);
1622 snew(newcgnr, at->nr);
1623 /* make index array to tell where the atoms go to when masses are inserted */
1624 snew(o2n, at->nr);
1625 for (i = 0; i < at->nr; i++)
1626 {
1627 o2n[i] = i;
1628 }
1629 /* make index to tell which residues were already processed */
1630 snew(bResProcessed, at->nres);
1631
1632 gmx_residuetype_init(&rt);
1633
1634 /* generate vsite constructions */
1635 /* loop over all atoms */
1636 resind = -1;
1637 for (i = 0; (i < at->nr); i++)
1638 {
1639 if (at->atom[i].resind != resind)
1640 {
1641 resind = at->atom[i].resind;
1642 resnm = *(at->resinfo[resind].name);
1643 }
1644 /* first check for aromatics to virtualize */
1645 /* don't waste our effort on DNA, water etc. */
1646 /* Only do the vsite aromatic stuff when we reach the
1647 * CA atom, since there might be an X2/X3 group on the
1648 * N-terminus that must be treated first.
1649 */
1650 if (bVsiteAromatics &&
1651 !strcmp(*(at->atomname[i]), "CA") &&
1652 !bResProcessed[resind] &&
1653 gmx_residuetype_is_protein(rt, *(at->resinfo[resind].name)) )
1654 {
1655 /* mark this residue */
1656 bResProcessed[resind] = TRUE;
1657 /* find out if this residue needs converting */
1658 whatres = NOTSET;
1659 for (j = 0; j < resNR && whatres == NOTSET; j++)
1660 {
1661
1662 cmplength = bPartial[j] ? strlen(resnm)-1 : strlen(resnm);
1663
1664 bFound = ((gmx_strncasecmp(resnm, resnms[j], cmplength) == 0) ||
1665 (gmx_strncasecmp(resnm, resnmsN[j], cmplength) == 0) ||
1666 (gmx_strncasecmp(resnm, resnmsC[j], cmplength) == 0));
1667
1668 if (bFound)
1669 {
1670 whatres = j;
1671 /* get atoms we will be needing for the conversion */
1672 nrfound = 0;
1673 for (k = 0; atnms[j][k]; k++)
1674 {
1675 ats[k] = NOTSET;
1676 for (m = i; m < at->nr && at->atom[m].resind == resind && ats[k] == NOTSET; m++)
1677 {
1678 if (gmx_strcasecmp(*(at->atomname[m]), atnms[j][k]) == 0)
1679 {
1680 ats[k] = m;
1681 nrfound++;
1682 }
1683 }
1684 }
1685
1686 /* now k is number of atom names in atnms[j] */
1687 if (j == resHIS)
1688 {
1689 needed = k-3;
1690 }
1691 else
1692 {
1693 needed = k;
1694 }
1695 if (nrfound < needed)
1696 {
1697 gmx_fatal(FARGS, "not enough atoms found (%d, need %d) in "
1698 "residue %s %d while\n "
1699 "generating aromatics virtual site construction",
1700 nrfound, needed, resnm, at->resinfo[resind].nr);
1701 }
1702 /* Advance overall atom counter */
1703 i++;
1704 }
1705 }
1706 /* the enums for every residue MUST correspond to atnms[residue] */
1707 switch (whatres)
1708 {
1709 case resPHE:
1710 if (debug)
1711 {
1712 fprintf(stderr, "PHE at %d\n", o2n[ats[0]]+1);
1713 }
1714 nvsite += gen_vsites_phe(at, vsite_type, plist, nrfound, ats, vsitetop, nvsitetop);
1715 break;
1716 case resTRP:
1717 if (debug)
1718 {
1719 fprintf(stderr, "TRP at %d\n", o2n[ats[0]]+1);
1720 }
1721 nvsite += gen_vsites_trp(atype, &newx, &newatom, &newatomname, &o2n,
1722 &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
1723 at, vsite_type, plist, nrfound, ats, add_shift, vsitetop, nvsitetop);
1724 break;
1725 case resTYR:
1726 if (debug)
1727 {
1728 fprintf(stderr, "TYR at %d\n", o2n[ats[0]]+1);
1729 }
1730 nvsite += gen_vsites_tyr(atype, &newx, &newatom, &newatomname, &o2n,
1731 &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
1732 at, vsite_type, plist, nrfound, ats, add_shift, vsitetop, nvsitetop);
1733 break;
1734 case resHIS:
1735 if (debug)
1736 {
1737 fprintf(stderr, "HIS at %d\n", o2n[ats[0]]+1);
1738 }
1739 nvsite += gen_vsites_his(at, vsite_type, plist, nrfound, ats, vsitetop, nvsitetop);
1740 break;
1741 case NOTSET:
1742 /* this means this residue won't be processed */
1743 break;
1744 default:
1745 gmx_fatal(FARGS, "DEATH HORROR in do_vsites (%s:%d)",
1746 __FILE__, __LINE__);
1747 } /* switch whatres */
1748 /* skip back to beginning of residue */
1749 while (i > 0 && at->atom[i-1].resind == resind)
1750 {
1751 i--;
1752 }
1753 } /* if bVsiteAromatics & is protein */
1754
1755 /* now process the rest of the hydrogens */
1756 /* only process hydrogen atoms which are not already set */
1757 if ( ((*vsite_type)[i] == NOTSET) && is_hydrogen(*(at->atomname[i])))
1758 {
1759 /* find heavy atom, count #bonds from it and #H atoms bound to it
1760 and return H atom numbers (Hatoms) and heavy atom numbers (heavies) */
1761 count_bonds(i, &plist[F_BONDS], at->atomname,
1762 &nrbonds, &nrHatoms, Hatoms, &Heavy, &nrheavies, heavies);
1763 /* get Heavy atom type */
1764 tpHeavy = get_atype(Heavy, at, nrtp, rtp, rt);
1765 strcpy(tpname, get_atomtype_name(tpHeavy, atype));
1766
1767 bWARNING = FALSE;
1768 bAddVsiteParam = TRUE;
1769 /* nested if's which check nrHatoms, nrbonds and atomname */
1770 if (nrHatoms == 1)
1771 {
1772 switch (nrbonds)
1773 {
1774 case 2: /* -O-H */
1775 (*vsite_type)[i] = F_BONDS;
1776 break;
1777 case 3: /* =CH-, -NH- or =NH+- */
1778 (*vsite_type)[i] = F_VSITE3FD;
1779 break;
1780 case 4: /* --CH- (tert) */
1781 /* The old type 4FD had stability issues, so
1782 * all new constructs should use 4FDN
1783 */
1784 (*vsite_type)[i] = F_VSITE4FDN;
1785
1786 /* Check parity of heavy atoms from coordinates */
1787 ai = Heavy;
1788 aj = heavies[0];
1789 ak = heavies[1];
1790 al = heavies[2];
1791 rvec_sub((*x)[aj], (*x)[ai], tmpmat[0]);
1792 rvec_sub((*x)[ak], (*x)[ai], tmpmat[1]);
1793 rvec_sub((*x)[al], (*x)[ai], tmpmat[2]);
1794
1795 if (det(tmpmat) > 0)
1796 {
1797 /* swap parity */
1798 heavies[1] = aj;
1799 heavies[0] = ak;
1800 }
1801
1802 break;
1803 default: /* nrbonds != 2, 3 or 4 */
1804 bWARNING = TRUE;
1805 }
1806
1807 }
1808 else if ( /*(nrHatoms == 2) && (nrbonds == 2) && REMOVED this test
1809 DvdS 19-01-04 */
1810 (gmx_strncasecmp(*at->atomname[Heavy], "OW", 2) == 0) )
1811 {
1812 bAddVsiteParam = FALSE; /* this is water: skip these hydrogens */
1813 if (bFirstWater)
1814 {
1815 bFirstWater = FALSE;
1816 if (debug)
1817 {
1818 fprintf(debug,
1819 "Not converting hydrogens in water to virtual sites\n");
1820 }
1821 }
1822 }
1823 else if ( (nrHatoms == 2) && (nrbonds == 4) )
1824 {
1825 /* -CH2- , -NH2+- */
1826 (*vsite_type)[Hatoms[0]] = F_VSITE3OUT;
1827 (*vsite_type)[Hatoms[1]] = -F_VSITE3OUT;
1828 }
1829 else
1830 {
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.
1833 */
1834 isN = planarN = FALSE;
1835 if ((nrHatoms == 2) && ((*at->atomname[Heavy])[0] == 'N'))
1836 {
1837 isN = TRUE;
1838 j = nitrogen_is_planar(vsiteconflist, nvsiteconf, tpname);
1839 if (j < 0)
1840 {
1841 gmx_fatal(FARGS, "No vsite database NH2 entry for type %s\n", tpname);
1842 }
1843 planarN = (j == 1);
1844 }
1845 if ( (nrHatoms == 2) && (nrbonds == 3) && ( !isN || planarN ) )
1846 {
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;
1850 }
1851 else if ( ( (nrHatoms == 2) && (nrbonds == 3) &&
1852 ( isN && !planarN ) ) ||
1853 ( (nrHatoms == 3) && (nrbonds == 4) ) )
1854 {
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 };
1858
1859 if (debug)
1860 {
1861 fprintf(stderr, "-XH3 or nonplanar NH2 group at %d\n", i+1);
1862 }
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++)
1867 {
1868 (*vsite_type)[Hatoms[j]] = Hat_vsite_type[j];
1869 }
1870 /* get dummy mass type from first char of heavy atom type (N or C) */
1871
1872 strcpy(nexttpname, get_atomtype_name(get_atype(heavies[0], at, nrtp, rtp, rt), atype));
1873 ch = get_dummymass_name(vsiteconflist, nvsiteconf, tpname, nexttpname);
1874
1875 if (ch == NULL)
1876 {
1877 if (ndb > 0)
1878 {
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);
1880 }
1881 else
1882 {
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);
1884 }
1885 }
1886 else
1887 {
1888 strcpy(name, ch);
1889 }
1890
1891 tpM = vsite_nm2type(name, atype);
1892 /* make space for 2 masses: shift all atoms starting with 'Heavy' */
1893 #define NMASS 2
1894 i0 = Heavy;
1895 ni0 = i0+nadd;
1896 if (debug)
1897 {
1898 fprintf(stderr, "Inserting %d dummy masses at %d\n", NMASS, o2n[i0]+1);
1899 }
1900 nadd += NMASS;
1901 for (j = i0; j < at->nr; j++)
1902 {
1903 o2n[j] = j+nadd;
1904 }
1905
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);
1911
1912 for (j = 0; j < NMASS; j++)
1913 {
1914 newatomname[at->nr+nadd-1-j] = NULL;
1915 }
1916
1917 /* calculate starting position for the masses */
1918 mHtot = 0;
1919 /* get atom masses, and set Heavy and Hatoms mass to zero */
1920 for (j = 0; j < nrHatoms; j++)
1921 {
1922 mHtot += get_amass(Hatoms[j], at, nrtp, rtp, rt);
1923 at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
1924 }
1925 mtot = mHtot + get_amass(Heavy, at, nrtp, rtp, rt);
1926 at->atom[Heavy].m = at->atom[Heavy].mB = 0;
1927 if (mHmult != 1.0)
1928 {
1929 mHtot *= mHmult;
1930 }
1931 fact2 = mHtot/mtot;
1932 fact = std::sqrt(fact2);
1933 /* generate vectors parallel and perpendicular to rotational axis:
1934 * rpar = Heavy -> Hcom
1935 * rperp = Hcom -> H1 */
1936 clear_rvec(rpar);
1937 for (j = 0; j < nrHatoms; j++)
1938 {
1939 rvec_inc(rpar, (*x)[Hatoms[j]]);
1940 }
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 */
1948 {
1949 rvec_add((*x)[Heavy], temp, newx[ni0+j]);
1950 }
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++)
1956 {
1957 /* make name: "M??#" or "M?#" (? is atomname, # is number) */
1958 name[0] = 'M';
1959 for (k = 0; (*at->atomname[Heavy])[k] && ( k < NMASS ); k++)
1960 {
1961 name[k+1] = (*at->atomname[Heavy])[k];
1962 }
1963 name[k+1] = atomnamesuffix[j];
1964 name[k+2] = '\0';
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];
1975 }
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);
1981
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,
1986 FALSE);
1987 for (j = 0; j < nrHatoms; j++)
1988 {
1989 add_vsite3_atoms(&plist[(*vsite_type)[Hatoms[j]]],
1990 Hatoms[j], heavies[0], add_shift+ni0, add_shift+ni0+1,
1991 Hat_SwapParity[j]);
1992 }
1993 #undef NMASS
1994 }
1995 else
1996 {
1997 bWARNING = TRUE;
1998 }
1999
2000 }
2001 if (bWARNING)
2002 {
2003 fprintf(stderr,
2004 "Warning: cannot convert atom %d %s (bound to a heavy atom "
2005 "%s with \n"
2006 " %d bonds and %d bound hydrogens atoms) to virtual site\n",
2007 i+1, *(at->atomname[i]), tpname, nrbonds, nrHatoms);
2008 }
2009 if (bAddVsiteParam)
2010 {
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++)
2017 {
2018 if (is_vsite((*vsite_type)[Hatoms[j]]))
2019 {
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;
2023 }
2024 }
2025 }
2026 nvsite += nrHatoms;
2027 if (debug)
2028 {
2029 fprintf(debug, "atom %d: ", o2n[i]+1);
2030 print_bonds(debug, o2n, nrHatoms, Hatoms, Heavy, nrheavies, heavies);
2031 }
2032 } /* if vsite NOTSET & is hydrogen */
2033
2034 } /* for i < at->nr */
2035
2036 gmx_residuetype_destroy(rt);
2037
2038 if (debug)
2039 {
2040 fprintf(debug, "Before inserting new atoms:\n");
2041 for (i = 0; i < at->nr; i++)
2042 {
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)",
2048 (*cgnr)[i],
2049 ((*vsite_type)[i] == NOTSET) ?
2050 "NOTSET" : interaction_function[(*vsite_type)[i]].name);
2051 }
2052 fprintf(debug, "new atoms to be inserted:\n");
2053 for (i = 0; i < at->nr+nadd; i++)
2054 {
2055 if (newatomname[i])
2056 {
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);
2062 }
2063 }
2064 }
2065
2066 /* add all original atoms to the new arrays, using o2n index array */
2067 for (i = 0; i < at->nr; i++)
2068 {
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]]);
2074 }
2075 /* throw away old atoms */
2076 sfree(at->atom);
2077 sfree(at->atomname);
2078 sfree(*vsite_type);
2079 sfree(*cgnr);
2080 sfree(*x);
2081 /* put in the new ones */
2082 at->nr += nadd;
2083 at->atom = newatom;
2084 at->atomname = newatomname;
2085 *vsite_type = newvsite_type;
2086 *cgnr = newcgnr;
2087 *x = newx;
2088 if (at->nr > add_shift)
2089 {
2090 gmx_fatal(FARGS, "Added impossible amount of dummy masses "
2091 "(%d on a total of %d atoms)\n", nadd, at->nr-nadd);
2092 }
2093
2094 if (debug)
2095 {
2096 fprintf(debug, "After inserting new atoms:\n");
2097 for (i = 0; i < at->nr; i++)
2098 {
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)",
2104 (*cgnr)[i],
2105 ((*vsite_type)[i] == NOTSET) ?
2106 "NOTSET" : interaction_function[(*vsite_type)[i]].name);
2107 }
2108 }
2109
2110 /* now renumber all the interactions because of the added atoms */
2111 for (ftype = 0; ftype < F_NRE; ftype++)
2112 {
2113 params = &(plist[ftype]);
2114 if (debug)
2115 {
2116 fprintf(debug, "Renumbering %d %s\n", params->nr,
2117 interaction_function[ftype].longname);
2118 }
2119 for (i = 0; i < params->nr; i++)
2120 {
2121 for (j = 0; j < NRAL(ftype); j++)
2122 {
2123 if (params->param[i].a[j] >= add_shift)
2124 {
2125 if (debug)
2126 {
2127 fprintf(debug, " [%d -> %d]", params->param[i].a[j],
2128 params->param[i].a[j]-add_shift);
2129 }
2130 params->param[i].a[j] = params->param[i].a[j]-add_shift;
2131 }
2132 else
2133 {
2134 if (debug)
2135 {
2136 fprintf(debug, " [%d -> %d]", params->param[i].a[j],
2137 o2n[params->param[i].a[j]]);
2138 }
2139 params->param[i].a[j] = o2n[params->param[i].a[j]];
2140 }
2141 }
2142 if (debug)
2143 {
2144 fprintf(debug, "\n");
2145 }
2146 }
2147 }
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++)
2151 {
2152 if (params->param[i].ai() > params->param[i].aj())
2153 {
2154 j = params->param[i].aj();
2155 params->param[i].aj() = params->param[i].ai();
2156 params->param[i].ai() = j;
2157 }
2158 }
2159
2160 /* clean up */
2161 sfree(o2n);
2162
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);
2167 }
2168
2169 void do_h_mass(t_params *psb, int vsite_type[], t_atoms *at, real mHmult,
2170 gmx_bool bDeuterate)
2171 {
2172 int i, j, a;
2173
2174 /* loop over all atoms */
2175 for (i = 0; i < at->nr; i++)
2176 {
2177 /* adjust masses if i is hydrogen and not a virtual site */
2178 if (!is_vsite(vsite_type[i]) && is_hydrogen(*(at->atomname[i])) )
2179 {
2180 /* find bonded heavy atom */
2181 a = NOTSET;
2182 for (j = 0; (j < psb->nr) && (a == NOTSET); j++)
2183 {
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()]) )
2187 {
2188 a = psb->param[j].aj();
2189 }
2190 else if ( (psb->param[j].aj() == i) &&
2191 !is_vsite(vsite_type[psb->param[j].ai()]) )
2192 {
2193 a = psb->param[j].ai();
2194 }
2195 }
2196 if (a == NOTSET)
2197 {
2198 gmx_fatal(FARGS, "Unbound hydrogen atom (%d) found while adjusting mass",
2199 i+1);
2200 }
2201
2202 /* adjust mass of i (hydrogen) with mHmult
2203 and correct mass of a (bonded atom) with same amount */
2204 if (!bDeuterate)
2205 {
2206 at->atom[a].m -= (mHmult-1.0)*at->atom[i].m;
2207 at->atom[a].mB -= (mHmult-1.0)*at->atom[i].m;
2208 }
2209 at->atom[i].m *= mHmult;
2210 at->atom[i].mB *= mHmult;
2211 }
2212 }
2213 }
The ultimate molecular dynamics simulation package