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