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