?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge branch 'master' of git://git.gromacs.org/gromacs
[gromacs/adressmacs.git] / src / kernel / gen_vsite.c
blobb3ce9f917cfc4a748ce0dd1fc61a9f69fefe968c
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 gmx_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(!gmx_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(!gmx_strcasecmp(dirstr,"HID"))
179 sprintf(dirstr,"HISA");
180 else if(!gmx_strcasecmp(dirstr,"HIE"))
181 sprintf(dirstr,"HISB");
182 else if(!gmx_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 && !gmx_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) && gmx_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 gmx_bool found=FALSE;
282 for(i=0;i<nvsiteconf && !found;i++) {
283 found=(!gmx_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 gmx_bool found=FALSE;
298 for(i=0;i<nvsiteconf && !found;i++) {
299 found=(!gmx_strcasecmp(vsiteconflist[i].atomtype,atom) &&
300 !gmx_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 && gmx_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 && gmx_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 gmx_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 = get_restp(*(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 gmx_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 = get_restp(*(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 gmx_bool bSwapParity;
487
488 for(i=0; i<nrHatoms; i++) {
489 ftype=vsite_type[Hatoms[i]];
490 /* Errors in setting the vsite_type should really be caugth earlier,
491 * because here it's not possible to print any useful error message.
492 * But it's still better to print a message than to segfault.
493 */
494 if (ftype == NOTSET) {
495 gmx_incons("Undetected error in setting up virtual sites");
496 }
497 bSwapParity = (ftype<0);
498 vsite_type[Hatoms[i]] = ftype = abs(ftype);
499 if (ftype == F_BONDS) {
500 if ( (nrheavies != 1) && (nrHatoms != 1) )
501 gmx_fatal(FARGS,"cannot make constraint in add_vsites for %d heavy "
502 "atoms and %d hydrogen atoms",nrheavies,nrHatoms);
503 my_add_param(&(plist[F_CONSTRNC]),Hatoms[i],heavies[0],NOTSET);
504 } else {
505 switch (ftype) {
506 case F_VSITE3:
507 case F_VSITE3FD:
508 case F_VSITE3OUT:
509 if (nrheavies < 2)
510 gmx_fatal(FARGS,"Not enough heavy atoms (%d) for %s (min 3)",
511 nrheavies+1,
512 interaction_function[vsite_type[Hatoms[i]]].name);
513 add_vsite3_atoms(&plist[ftype],Hatoms[i],Heavy,heavies[0],heavies[1],
514 bSwapParity);
515 break;
516 case F_VSITE3FAD: {
517 if (nrheavies > 1)
518 moreheavy=heavies[1];
519 else {
520 /* find more heavy atoms */
521 other=moreheavy=NOTSET;
522 for(j=0; (j<plist[F_BONDS].nr) && (moreheavy==NOTSET); j++) {
523 if (plist[F_BONDS].param[j].AI==heavies[0])
524 other=plist[F_BONDS].param[j].AJ;
525 else if (plist[F_BONDS].param[j].AJ==heavies[0])
526 other=plist[F_BONDS].param[j].AI;
527 if ( (other != NOTSET) && (other != Heavy) )
528 moreheavy=other;
529 }
530 if (moreheavy==NOTSET)
531 gmx_fatal(FARGS,"Unbound molecule part %d-%d",Heavy+1,Hatoms[0]+1);
532 }
533 add_vsite3_atoms(&plist[ftype],Hatoms[i],Heavy,heavies[0],moreheavy,
534 bSwapParity);
535 break;
536 }
537 case F_VSITE4FD:
538 case F_VSITE4FDN:
539 if (nrheavies < 3)
540 {
541 gmx_fatal(FARGS,"Not enough heavy atoms (%d) for %s (min 4)",
542 nrheavies+1,
543 interaction_function[vsite_type[Hatoms[i]]].name);
544 }
545 add_vsite4_atoms(&plist[ftype],
546 Hatoms[0], Heavy, heavies[0], heavies[1], heavies[2]);
547 break;
548
549 default:
550 gmx_fatal(FARGS,"can't use add_vsites for interaction function %s",
551 interaction_function[vsite_type[Hatoms[i]]].name);
552 } /* switch ftype */
553 } /* else */
554 } /* for i */
555 }
556
557 #define ANGLE_6RING (DEG2RAD*120)
558
559 /* cosine rule: a^2 = b^2 + c^2 - 2 b c cos(alpha) */
560 /* get a^2 when a, b and alpha are given: */
561 #define cosrule(b,c,alpha) ( sqr(b) + sqr(c) - 2*b*c*cos(alpha) )
562 /* get cos(alpha) when a, b and c are given: */
563 #define acosrule(a,b,c) ( (sqr(b)+sqr(c)-sqr(a))/(2*b*c) )
564
565 static int gen_vsites_6ring(t_atoms *at, int *vsite_type[], t_params plist[],
566 int nrfound, int *ats, real bond_cc, real bond_ch,
567 real xcom, real ycom, gmx_bool bDoZ)
568 {
569 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
570 enum { atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
571 atCZ, atHZ, atNR };
572
573 int i,nvsite;
574 real a,b,dCGCE,tmp1,tmp2,mtot,mG,mrest;
575 real xCG,yCG,xCE1,yCE1,xCE2,yCE2;
576 /* CG, CE1 and CE2 stay and each get a part of the total mass,
577 * so the c-o-m stays the same.
578 */
579
580 if (bDoZ) {
581 if (atNR != nrfound)
582 gmx_incons("Generating vsites on 6-rings");
583 }
584
585 /* constraints between CG, CE1 and CE2: */
586 dCGCE = sqrt( cosrule(bond_cc,bond_cc,ANGLE_6RING) );
587 my_add_param(&(plist[F_CONSTRNC]),ats[atCG] ,ats[atCE1],dCGCE);
588 my_add_param(&(plist[F_CONSTRNC]),ats[atCG] ,ats[atCE2],dCGCE);
589 my_add_param(&(plist[F_CONSTRNC]),ats[atCE1],ats[atCE2],dCGCE);
590
591 /* rest will be vsite3 */
592 mtot=0;
593 nvsite=0;
594 for(i=0; i< (bDoZ ? atNR : atHZ) ; i++) {
595 mtot+=at->atom[ats[i]].m;
596 if ( i!=atCG && i!=atCE1 && i!=atCE2 && (bDoZ || (i!=atHZ && i!=atCZ) ) ) {
597 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
598 (*vsite_type)[ats[i]]=F_VSITE3;
599 nvsite++;
600 }
601 }
602 /* Distribute mass so center-of-mass stays the same.
603 * The center-of-mass in the call is defined with x=0 at
604 * the CE1-CE2 bond and y=0 at the line from CG to the middle of CE1-CE2 bond.
605 */
606 xCG=-bond_cc+bond_cc*cos(ANGLE_6RING);
607 yCG=0;
608 xCE1=0;
609 yCE1=bond_cc*sin(0.5*ANGLE_6RING);
610 xCE2=0;
611 yCE2=-bond_cc*sin(0.5*ANGLE_6RING);
612
613 mG = at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = xcom*mtot/xCG;
614 mrest=mtot-mG;
615 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB =
616 at->atom[ats[atCE2]].m = at->atom[ats[atCE2]].mB = mrest / 2;
617
618 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
619 tmp1 = dCGCE*sin(ANGLE_6RING*0.5);
620 tmp2 = bond_cc*cos(0.5*ANGLE_6RING) + tmp1;
621 tmp1 *= 2;
622 a = b = - bond_ch / tmp1;
623 /* HE1 and HE2: */
624 add_vsite3_param(&plist[F_VSITE3],
625 ats[atHE1],ats[atCE1],ats[atCE2],ats[atCG], a,b);
626 add_vsite3_param(&plist[F_VSITE3],
627 ats[atHE2],ats[atCE2],ats[atCE1],ats[atCG], a,b);
628 /* CD1, CD2 and CZ: */
629 a = b = tmp2 / tmp1;
630 add_vsite3_param(&plist[F_VSITE3],
631 ats[atCD1],ats[atCE2],ats[atCE1],ats[atCG], a,b);
632 add_vsite3_param(&plist[F_VSITE3],
633 ats[atCD2],ats[atCE1],ats[atCE2],ats[atCG], a,b);
634 if (bDoZ)
635 add_vsite3_param(&plist[F_VSITE3],
636 ats[atCZ], ats[atCG], ats[atCE1],ats[atCE2],a,b);
637 /* HD1, HD2 and HZ: */
638 a = b = ( bond_ch + tmp2 ) / tmp1;
639 add_vsite3_param(&plist[F_VSITE3],
640 ats[atHD1],ats[atCE2],ats[atCE1],ats[atCG], a,b);
641 add_vsite3_param(&plist[F_VSITE3],
642 ats[atHD2],ats[atCE1],ats[atCE2],ats[atCG], a,b);
643 if (bDoZ)
644 add_vsite3_param(&plist[F_VSITE3],
645 ats[atHZ], ats[atCG], ats[atCE1],ats[atCE2],a,b);
646
647 return nvsite;
648 }
649
650 static int gen_vsites_phe(t_atoms *at, int *vsite_type[], t_params plist[],
651 int nrfound, int *ats, t_vsitetop *vsitetop, int nvsitetop)
652 {
653 real bond_cc,bond_ch;
654 real xcom,ycom,mtot;
655 int i;
656 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
657 enum { atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
658 atCZ, atHZ, atNR };
659 real x[atNR],y[atNR];
660 /* Aromatic rings have 6-fold symmetry, so we only need one bond length.
661 * (angle is always 120 degrees).
662 */
663 bond_cc=get_ddb_bond(vsitetop,nvsitetop,"PHE","CD1","CE1");
664 bond_ch=get_ddb_bond(vsitetop,nvsitetop,"PHE","CD1","HD1");
665
666 x[atCG]=-bond_cc+bond_cc*cos(ANGLE_6RING);
667 y[atCG]=0;
668 x[atCD1]=-bond_cc;
669 y[atCD1]=bond_cc*sin(0.5*ANGLE_6RING);
670 x[atHD1]=x[atCD1]+bond_ch*cos(ANGLE_6RING);
671 y[atHD1]=y[atCD1]+bond_ch*sin(ANGLE_6RING);
672 x[atCE1]=0;
673 y[atCE1]=y[atCD1];
674 x[atHE1]=x[atCE1]-bond_ch*cos(ANGLE_6RING);
675 y[atHE1]=y[atCE1]+bond_ch*sin(ANGLE_6RING);
676 x[atCD2]=x[atCD1];
677 y[atCD2]=-y[atCD1];
678 x[atHD2]=x[atHD1];
679 y[atHD2]=-y[atHD1];
680 x[atCE2]=x[atCE1];
681 y[atCE2]=-y[atCE1];
682 x[atHE2]=x[atHE1];
683 y[atHE2]=-y[atHE1];
684 x[atCZ]=bond_cc*cos(0.5*ANGLE_6RING);
685 y[atCZ]=0;
686 x[atHZ]=x[atCZ]+bond_ch;
687 y[atHZ]=0;
688
689 xcom=ycom=mtot=0;
690 for(i=0;i<atNR;i++) {
691 xcom+=x[i]*at->atom[ats[i]].m;
692 ycom+=y[i]*at->atom[ats[i]].m;
693 mtot+=at->atom[ats[i]].m;
694 }
695 xcom/=mtot;
696 ycom/=mtot;
697
698 return gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, ycom, TRUE);
699 }
700
701 static void calc_vsite3_param(real xd,real yd,real xi,real yi,real xj,real yj,
702 real xk, real yk, real *a, real *b)
703 {
704 /* determine parameters by solving the equation system, since we know the
705 * virtual site coordinates here.
706 */
707 real dx_ij,dx_ik,dy_ij,dy_ik;
708 real b_ij,b_ik;
709
710 dx_ij=xj-xi;
711 dy_ij=yj-yi;
712 dx_ik=xk-xi;
713 dy_ik=yk-yi;
714 b_ij=sqrt(dx_ij*dx_ij+dy_ij*dy_ij);
715 b_ik=sqrt(dx_ik*dx_ik+dy_ik*dy_ik);
716
717 *a = ( (xd-xi)*dy_ik - dx_ik*(yd-yi) ) / (dx_ij*dy_ik - dx_ik*dy_ij);
718 *b = ( yd - yi - (*a)*dy_ij ) / dy_ik;
719 }
720
721
722 static int gen_vsites_trp(gpp_atomtype_t atype, rvec *newx[],
723 t_atom *newatom[], char ***newatomname[],
724 int *o2n[], int *newvsite_type[], int *newcgnr[],
725 t_symtab *symtab, int *nadd, rvec x[], int *cgnr[],
726 t_atoms *at, int *vsite_type[], t_params plist[],
727 int nrfound, int *ats, int add_shift,
728 t_vsitetop *vsitetop, int nvsitetop)
729 {
730 #define NMASS 2
731 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
732 enum {
733 atCB, atCG, atCD1, atHD1, atCD2, atNE1, atHE1, atCE2, atCE3, atHE3,
734 atCZ2, atHZ2, atCZ3, atHZ3, atCH2, atHH2, atNR };
735 /* weights for determining the COM's of both rings (M1 and M2): */
736 real mw[NMASS][atNR] = {
737 { 0, 1, 1, 1, 0.5, 1, 1, 0.5, 0, 0,
738 0, 0, 0, 0, 0, 0 },
739 { 0, 0, 0, 0, 0.5, 0, 0, 0.5, 1, 1,
740 1, 1, 1, 1, 1, 1 }
741 };
742
743 real xi[atNR],yi[atNR];
744 real xcom[NMASS],ycom[NMASS],I,alpha;
745 real lineA,lineB,dist;
746 real b_CD2_CE2,b_NE1_CE2,b_CG_CD2,b_CH2_HH2,b_CE2_CZ2;
747 real b_NE1_HE1,b_CD2_CE3,b_CE3_CZ3,b_CB_CG;
748 real b_CZ2_CH2,b_CZ2_HZ2,b_CD1_HD1,b_CE3_HE3;
749 real b_CG_CD1,b_CZ3_HZ3;
750 real a_NE1_CE2_CD2,a_CE2_CD2_CG,a_CB_CG_CD2,a_CE2_CD2_CE3;
751 real a_CB_CG_CD1,a_CD2_CG_CD1,a_CE2_CZ2_HZ2,a_CZ2_CH2_HH2;
752 real a_CD2_CE2_CZ2,a_CD2_CE3_CZ3,a_CE3_CZ3_HZ3,a_CG_CD1_HD1;
753 real a_CE2_CZ2_CH2,a_HE1_NE1_CE2,a_CD2_CE3_HE3;
754 real xM[NMASS];
755 int atM[NMASS],tpM,i,i0,j,nvsite;
756 real mwtot,mtot,mM[NMASS],dCBM1,dCBM2,dM1M2;
757 real a,b,c[MAXFORCEPARAM];
758 rvec r_ij,r_ik,t1,t2;
759 char name[10];
760
761 if (atNR != nrfound)
762 gmx_incons("atom types in gen_vsites_trp");
763 /* Get geometry from database */
764 b_CD2_CE2=get_ddb_bond(vsitetop,nvsitetop,"TRP","CD2","CE2");
765 b_NE1_CE2=get_ddb_bond(vsitetop,nvsitetop,"TRP","NE1","CE2");
766 b_CG_CD1=get_ddb_bond(vsitetop,nvsitetop,"TRP","CG","CD1");
767 b_CG_CD2=get_ddb_bond(vsitetop,nvsitetop,"TRP","CG","CD2");
768 b_CB_CG=get_ddb_bond(vsitetop,nvsitetop,"TRP","CB","CG");
769 b_CE2_CZ2=get_ddb_bond(vsitetop,nvsitetop,"TRP","CE2","CZ2");
770 b_CD2_CE3=get_ddb_bond(vsitetop,nvsitetop,"TRP","CD2","CE3");
771 b_CE3_CZ3=get_ddb_bond(vsitetop,nvsitetop,"TRP","CE3","CZ3");
772 b_CZ2_CH2=get_ddb_bond(vsitetop,nvsitetop,"TRP","CZ2","CH2");
773
774 b_CD1_HD1=get_ddb_bond(vsitetop,nvsitetop,"TRP","CD1","HD1");
775 b_CZ2_HZ2=get_ddb_bond(vsitetop,nvsitetop,"TRP","CZ2","HZ2");
776 b_NE1_HE1=get_ddb_bond(vsitetop,nvsitetop,"TRP","NE1","HE1");
777 b_CH2_HH2=get_ddb_bond(vsitetop,nvsitetop,"TRP","CH2","HH2");
778 b_CE3_HE3=get_ddb_bond(vsitetop,nvsitetop,"TRP","CE3","HE3");
779 b_CZ3_HZ3=get_ddb_bond(vsitetop,nvsitetop,"TRP","CZ3","HZ3");
780
781 a_NE1_CE2_CD2=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","NE1","CE2","CD2");
782 a_CE2_CD2_CG=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","CE2","CD2","CG");
783 a_CB_CG_CD2=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","CB","CG","CD2");
784 a_CD2_CG_CD1=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","CD2","CG","CD1");
785 a_CB_CG_CD1=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","CB","CG","CD1");
786
787 a_CE2_CD2_CE3=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","CE2","CD2","CE3");
788 a_CD2_CE2_CZ2=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","CD2","CE2","CZ2");
789 a_CD2_CE3_CZ3=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","CD2","CE3","CZ3");
790 a_CE3_CZ3_HZ3=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","CE3","CZ3","HZ3");
791 a_CZ2_CH2_HH2=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","CZ2","CH2","HH2");
792 a_CE2_CZ2_HZ2=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","CE2","CZ2","HZ2");
793 a_CE2_CZ2_CH2=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","CE2","CZ2","CH2");
794 a_CG_CD1_HD1=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","CG","CD1","HD1");
795 a_HE1_NE1_CE2=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","HE1","NE1","CE2");
796 a_CD2_CE3_HE3=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TRP","CD2","CE3","HE3");
797
798 /* Calculate local coordinates.
799 * y-axis (x=0) is the bond CD2-CE2.
800 * x-axis (y=0) is perpendicular to the bond CD2-CE2 and
801 * intersects the middle of the bond.
802 */
803 xi[atCD2]=0;
804 yi[atCD2]=-0.5*b_CD2_CE2;
805
806 xi[atCE2]=0;
807 yi[atCE2]=0.5*b_CD2_CE2;
808
809 xi[atNE1]=-b_NE1_CE2*sin(a_NE1_CE2_CD2);
810 yi[atNE1]=yi[atCE2]-b_NE1_CE2*cos(a_NE1_CE2_CD2);
811
812 xi[atCG]=-b_CG_CD2*sin(a_CE2_CD2_CG);
813 yi[atCG]=yi[atCD2]+b_CG_CD2*cos(a_CE2_CD2_CG);
814
815 alpha = a_CE2_CD2_CG + M_PI - a_CB_CG_CD2;
816 xi[atCB]=xi[atCG]-b_CB_CG*sin(alpha);
817 yi[atCB]=yi[atCG]+b_CB_CG*cos(alpha);
818
819 alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - M_PI;
820 xi[atCD1]=xi[atCG]-b_CG_CD1*sin(alpha);
821 yi[atCD1]=yi[atCG]+b_CG_CD1*cos(alpha);
822
823 xi[atCE3]=b_CD2_CE3*sin(a_CE2_CD2_CE3);
824 yi[atCE3]=yi[atCD2]+b_CD2_CE3*cos(a_CE2_CD2_CE3);
825
826 xi[atCZ2]=b_CE2_CZ2*sin(a_CD2_CE2_CZ2);
827 yi[atCZ2]=yi[atCE2]-b_CE2_CZ2*cos(a_CD2_CE2_CZ2);
828
829 alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - M_PI;
830 xi[atCZ3]=xi[atCE3]+b_CE3_CZ3*sin(alpha);
831 yi[atCZ3]=yi[atCE3]+b_CE3_CZ3*cos(alpha);
832
833 alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - M_PI;
834 xi[atCH2]=xi[atCZ2]+b_CZ2_CH2*sin(alpha);
835 yi[atCH2]=yi[atCZ2]-b_CZ2_CH2*cos(alpha);
836
837 /* hydrogens */
838 alpha = a_CE2_CD2_CG + a_CD2_CG_CD1 - a_CG_CD1_HD1;
839 xi[atHD1]=xi[atCD1]-b_CD1_HD1*sin(alpha);
840 yi[atHD1]=yi[atCD1]+b_CD1_HD1*cos(alpha);
841
842 alpha = a_NE1_CE2_CD2 + M_PI - a_HE1_NE1_CE2;
843 xi[atHE1]=xi[atNE1]-b_NE1_HE1*sin(alpha);
844 yi[atHE1]=yi[atNE1]-b_NE1_HE1*cos(alpha);
845
846 alpha = a_CE2_CD2_CE3 + M_PI - a_CD2_CE3_HE3;
847 xi[atHE3]=xi[atCE3]+b_CE3_HE3*sin(alpha);
848 yi[atHE3]=yi[atCE3]+b_CE3_HE3*cos(alpha);
849
850 alpha = a_CD2_CE2_CZ2 + M_PI - a_CE2_CZ2_HZ2;
851 xi[atHZ2]=xi[atCZ2]+b_CZ2_HZ2*sin(alpha);
852 yi[atHZ2]=yi[atCZ2]-b_CZ2_HZ2*cos(alpha);
853
854 alpha = a_CD2_CE2_CZ2 + a_CE2_CZ2_CH2 - a_CZ2_CH2_HH2;
855 xi[atHZ3]=xi[atCZ3]+b_CZ3_HZ3*sin(alpha);
856 yi[atHZ3]=yi[atCZ3]+b_CZ3_HZ3*cos(alpha);
857
858 alpha = a_CE2_CD2_CE3 + a_CD2_CE3_CZ3 - a_CE3_CZ3_HZ3;
859 xi[atHH2]=xi[atCH2]+b_CH2_HH2*sin(alpha);
860 yi[atHH2]=yi[atCH2]-b_CH2_HH2*cos(alpha);
861
862 /* Determine coeff. for the line CB-CG */
863 lineA=(yi[atCB]-yi[atCG])/(xi[atCB]-xi[atCG]);
864 lineB=yi[atCG]-lineA*xi[atCG];
865
866 /* Calculate masses for each ring and put it on the dummy masses */
867 for(j=0; j<NMASS; j++)
868 mM[j]=xcom[j]=ycom[j]=0;
869 for(i=0; i<atNR; i++) {
870 if (i!=atCB) {
871 for(j=0; j<NMASS; j++) {
872 mM[j] += mw[j][i] * at->atom[ats[i]].m;
873 xcom[j] += xi[i] * mw[j][i] * at->atom[ats[i]].m;
874 ycom[j] += yi[i] * mw[j][i] * at->atom[ats[i]].m;
875 }
876 }
877 }
878 for(j=0; j<NMASS; j++) {
879 xcom[j]/=mM[j];
880 ycom[j]/=mM[j];
881 }
882
883 /* get dummy mass type */
884 tpM=vsite_nm2type("MW",atype);
885 /* make space for 2 masses: shift all atoms starting with CB */
886 i0=ats[atCB];
887 for(j=0; j<NMASS; j++)
888 atM[j]=i0+*nadd+j;
889 if (debug)
890 fprintf(stderr,"Inserting %d dummy masses at %d\n",NMASS,(*o2n)[i0]+1);
891 *nadd+=NMASS;
892 for(j=i0; j<at->nr; j++)
893 (*o2n)[j]=j+*nadd;
894 srenew(*newx,at->nr+*nadd);
895 srenew(*newatom,at->nr+*nadd);
896 srenew(*newatomname,at->nr+*nadd);
897 srenew(*newvsite_type,at->nr+*nadd);
898 srenew(*newcgnr,at->nr+*nadd);
899 for(j=0; j<NMASS; j++)
900 (*newatomname)[at->nr+*nadd-1-j]=NULL;
901
902 /* Dummy masses will be placed at the center-of-mass in each ring. */
903
904 /* calc initial position for dummy masses in real (non-local) coordinates.
905 * Cheat by using the routine to calculate virtual site parameters. It is
906 * much easier when we have the coordinates expressed in terms of
907 * CB, CG, CD2.
908 */
909 rvec_sub(x[ats[atCB]],x[ats[atCG]],r_ij);
910 rvec_sub(x[ats[atCD2]],x[ats[atCG]],r_ik);
911 calc_vsite3_param(xcom[0],ycom[0],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[0]]);
917
918 calc_vsite3_param(xcom[1],ycom[1],xi[atCG],yi[atCG],xi[atCB],yi[atCB],
919 xi[atCD2],yi[atCD2],&a,&b);
920 svmul(a,r_ij,t1);
921 svmul(b,r_ik,t2);
922 rvec_add(t1,t2,t1);
923 rvec_add(t1,x[ats[atCG]],(*newx)[atM[1]]);
924
925 /* set parameters for the masses */
926 for(j=0; j<NMASS; j++) {
927 sprintf(name,"MW%d",j+1);
928 (*newatomname) [atM[j]] = put_symtab(symtab,name);
929 (*newatom) [atM[j]].m = (*newatom)[atM[j]].mB = mM[j];
930 (*newatom) [atM[j]].q = (*newatom)[atM[j]].qB = 0.0;
931 (*newatom) [atM[j]].type = (*newatom)[atM[j]].typeB = tpM;
932 (*newatom) [atM[j]].ptype = eptAtom;
933 (*newatom) [atM[j]].resind= at->atom[i0].resind;
934 (*newvsite_type)[atM[j]] = NOTSET;
935 (*newcgnr) [atM[j]] = (*cgnr)[i0];
936 }
937 /* renumber cgnr: */
938 for(i=i0; i<at->nr; i++)
939 (*cgnr)[i]++;
940
941 /* constraints between CB, M1 and M2 */
942 /* 'add_shift' says which atoms won't be renumbered afterwards */
943 dCBM1 = sqrt( sqr(xcom[0]-xi[atCB]) + sqr(ycom[0]-yi[atCB]) );
944 dM1M2 = sqrt( sqr(xcom[0]-xcom[1]) + sqr(ycom[0]-ycom[1]) );
945 dCBM2 = sqrt( sqr(xcom[1]-xi[atCB]) + sqr(ycom[1]-yi[atCB]) );
946 my_add_param(&(plist[F_CONSTRNC]),ats[atCB], add_shift+atM[0],dCBM1);
947 my_add_param(&(plist[F_CONSTRNC]),ats[atCB], add_shift+atM[1],dCBM2);
948 my_add_param(&(plist[F_CONSTRNC]),add_shift+atM[0],add_shift+atM[1],dM1M2);
949
950 /* rest will be vsite3 */
951 nvsite=0;
952 for(i=0; i<atNR; i++)
953 if (i!=atCB) {
954 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
955 (*vsite_type)[ats[i]] = F_VSITE3;
956 nvsite++;
957 }
958
959 /* now define all vsites from M1, M2, CB, ie:
960 r_d = r_M1 + a r_M1_M2 + b r_M1_CB */
961 for(i=0; i<atNR; i++)
962 if ( (*vsite_type)[ats[i]] == F_VSITE3 ) {
963 calc_vsite3_param(xi[i],yi[i],xcom[0],ycom[0],xcom[1],ycom[1],xi[atCB],yi[atCB],&a,&b);
964 add_vsite3_param(&plist[F_VSITE3],
965 ats[i],add_shift+atM[0],add_shift+atM[1],ats[atCB],a,b);
966 }
967 return nvsite;
968 #undef NMASS
969 }
970
971
972 static int gen_vsites_tyr(gpp_atomtype_t atype, rvec *newx[],
973 t_atom *newatom[], char ***newatomname[],
974 int *o2n[], int *newvsite_type[], int *newcgnr[],
975 t_symtab *symtab, int *nadd, rvec x[], int *cgnr[],
976 t_atoms *at, int *vsite_type[], t_params plist[],
977 int nrfound, int *ats, int add_shift,
978 t_vsitetop *vsitetop, int nvsitetop)
979 {
980 int nvsite,i,i0,j,atM,tpM;
981 real dCGCE,dCEOH,dCGM,tmp1,a,b;
982 real bond_cc,bond_ch,bond_co,bond_oh,angle_coh;
983 real xcom,ycom,mtot;
984 real vmass,vdist,mM;
985 rvec r1;
986 char name[10];
987
988 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
989 enum { atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
990 atCZ, atOH, atHH, atNR };
991 real xi[atNR],yi[atNR];
992 /* CG, CE1, CE2 (as in general 6-ring) and OH and HH stay,
993 rest gets virtualized.
994 Now we have two linked triangles with one improper keeping them flat */
995 if (atNR != nrfound)
996 gmx_incons("Number of atom types in gen_vsites_tyr");
997
998 /* Aromatic rings have 6-fold symmetry, so we only need one bond length
999 * for the ring part (angle is always 120 degrees).
1000 */
1001 bond_cc=get_ddb_bond(vsitetop,nvsitetop,"TYR","CD1","CE1");
1002 bond_ch=get_ddb_bond(vsitetop,nvsitetop,"TYR","CD1","HD1");
1003 bond_co=get_ddb_bond(vsitetop,nvsitetop,"TYR","CZ","OH");
1004 bond_oh=get_ddb_bond(vsitetop,nvsitetop,"TYR","OH","HH");
1005 angle_coh=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,"TYR","CZ","OH","HH");
1006
1007 xi[atCG]=-bond_cc+bond_cc*cos(ANGLE_6RING);
1008 yi[atCG]=0;
1009 xi[atCD1]=-bond_cc;
1010 yi[atCD1]=bond_cc*sin(0.5*ANGLE_6RING);
1011 xi[atHD1]=xi[atCD1]+bond_ch*cos(ANGLE_6RING);
1012 yi[atHD1]=yi[atCD1]+bond_ch*sin(ANGLE_6RING);
1013 xi[atCE1]=0;
1014 yi[atCE1]=yi[atCD1];
1015 xi[atHE1]=xi[atCE1]-bond_ch*cos(ANGLE_6RING);
1016 yi[atHE1]=yi[atCE1]+bond_ch*sin(ANGLE_6RING);
1017 xi[atCD2]=xi[atCD1];
1018 yi[atCD2]=-yi[atCD1];
1019 xi[atHD2]=xi[atHD1];
1020 yi[atHD2]=-yi[atHD1];
1021 xi[atCE2]=xi[atCE1];
1022 yi[atCE2]=-yi[atCE1];
1023 xi[atHE2]=xi[atHE1];
1024 yi[atHE2]=-yi[atHE1];
1025 xi[atCZ]=bond_cc*cos(0.5*ANGLE_6RING);
1026 yi[atCZ]=0;
1027 xi[atOH]=xi[atCZ]+bond_co;
1028 yi[atOH]=0;
1029
1030 xcom=ycom=mtot=0;
1031 for(i=0;i<atOH;i++) {
1032 xcom+=xi[i]*at->atom[ats[i]].m;
1033 ycom+=yi[i]*at->atom[ats[i]].m;
1034 mtot+=at->atom[ats[i]].m;
1035 }
1036 xcom/=mtot;
1037 ycom/=mtot;
1038
1039 /* first do 6 ring as default,
1040 except CZ (we'll do that different) and HZ (we don't have that): */
1041 nvsite = gen_vsites_6ring(at, vsite_type, plist, nrfound, ats, bond_cc, bond_ch, xcom, ycom, FALSE);
1042
1043 /* then construct CZ from the 2nd triangle */
1044 /* vsite3 construction: r_d = r_i + a r_ij + b r_ik */
1045 a = b = 0.5 * bond_co / ( bond_co - bond_cc*cos(ANGLE_6RING) );
1046 add_vsite3_param(&plist[F_VSITE3],
1047 ats[atCZ],ats[atOH],ats[atCE1],ats[atCE2],a,b);
1048 at->atom[ats[atCZ]].m = at->atom[ats[atCZ]].mB = 0;
1049
1050 /* constraints between CE1, CE2 and OH */
1051 dCGCE = sqrt( cosrule(bond_cc,bond_cc,ANGLE_6RING) );
1052 dCEOH = sqrt( cosrule(bond_cc,bond_co,ANGLE_6RING) );
1053 my_add_param(&(plist[F_CONSTRNC]),ats[atCE1],ats[atOH],dCEOH);
1054 my_add_param(&(plist[F_CONSTRNC]),ats[atCE2],ats[atOH],dCEOH);
1055
1056 /* We also want to constrain the angle C-O-H, but since CZ is constructed
1057 * we need to introduce a constraint to CG.
1058 * CG is much further away, so that will lead to instabilities in LINCS
1059 * when we constrain both CG-HH and OH-HH distances. Instead of requiring
1060 * the use of lincs_order=8 we introduce a dummy mass three times further
1061 * away from OH than HH. The mass is accordingly a third, with the remaining
1062 * 2/3 moved to OH. This shouldnt cause any problems since the forces will
1063 * apply to the HH constructed atom and not directly on the virtual mass.
1064 */
1065
1066 vdist=2.0*bond_oh;
1067 mM=at->atom[ats[atHH]].m/2.0;
1068 at->atom[ats[atOH]].m+=mM; /* add 1/2 of original H mass */
1069 at->atom[ats[atOH]].mB+=mM; /* add 1/2 of original H mass */
1070 at->atom[ats[atHH]].m=at->atom[ats[atHH]].mB=0;
1071
1072 /* get dummy mass type */
1073 tpM=vsite_nm2type("MW",atype);
1074 /* make space for 1 mass: shift HH only */
1075 i0=ats[atHH];
1076 atM=i0+*nadd;
1077 if (debug)
1078 fprintf(stderr,"Inserting 1 dummy mass at %d\n",(*o2n)[i0]+1);
1079 (*nadd)++;
1080 for(j=i0; j<at->nr; j++)
1081 (*o2n)[j]=j+*nadd;
1082 srenew(*newx,at->nr+*nadd);
1083 srenew(*newatom,at->nr+*nadd);
1084 srenew(*newatomname,at->nr+*nadd);
1085 srenew(*newvsite_type,at->nr+*nadd);
1086 srenew(*newcgnr,at->nr+*nadd);
1087 (*newatomname)[at->nr+*nadd-1]=NULL;
1088
1089 /* Calc the dummy mass initial position */
1090 rvec_sub(x[ats[atHH]],x[ats[atOH]],r1);
1091 svmul(2.0,r1,r1);
1092 rvec_add(r1,x[ats[atHH]],(*newx)[atM]);
1093
1094 strcpy(name,"MW1");
1095 (*newatomname) [atM] = put_symtab(symtab,name);
1096 (*newatom) [atM].m = (*newatom)[atM].mB = mM;
1097 (*newatom) [atM].q = (*newatom)[atM].qB = 0.0;
1098 (*newatom) [atM].type = (*newatom)[atM].typeB = tpM;
1099 (*newatom) [atM].ptype = eptAtom;
1100 (*newatom) [atM].resind= at->atom[i0].resind;
1101 (*newvsite_type)[atM] = NOTSET;
1102 (*newcgnr) [atM] = (*cgnr)[i0];
1103 /* renumber cgnr: */
1104 for(i=i0; i<at->nr; i++)
1105 (*cgnr)[i]++;
1106
1107 (*vsite_type)[ats[atHH]] = F_VSITE2;
1108 nvsite++;
1109 /* assume we also want the COH angle constrained: */
1110 tmp1 = bond_cc*cos(0.5*ANGLE_6RING) + dCGCE*sin(ANGLE_6RING*0.5) + bond_co;
1111 dCGM = sqrt( cosrule(tmp1,vdist,angle_coh) );
1112 my_add_param(&(plist[F_CONSTRNC]),ats[atCG],add_shift+atM,dCGM);
1113 my_add_param(&(plist[F_CONSTRNC]),ats[atOH],add_shift+atM,vdist);
1114
1115 add_vsite2_param(&plist[F_VSITE2],
1116 ats[atHH],ats[atOH],add_shift+atM,1.0/2.0);
1117 return nvsite;
1118 }
1119
1120 static int gen_vsites_his(t_atoms *at, int *vsite_type[], t_params plist[],
1121 int nrfound, int *ats, t_vsitetop *vsitetop, int nvsitetop)
1122 {
1123 int nvsite,i;
1124 real a,b,alpha,dCGCE1,dCGNE2;
1125 real sinalpha,cosalpha;
1126 real xcom,ycom,mtot;
1127 real mG,mrest,mCE1,mNE2;
1128 real b_CG_ND1,b_ND1_CE1,b_CE1_NE2,b_CG_CD2,b_CD2_NE2;
1129 real b_ND1_HD1,b_NE2_HE2,b_CE1_HE1,b_CD2_HD2;
1130 real a_CG_ND1_CE1,a_CG_CD2_NE2,a_ND1_CE1_NE2,a_CE1_NE2_CD2;
1131 real a_NE2_CE1_HE1,a_NE2_CD2_HD2,a_CE1_ND1_HD1,a_CE1_NE2_HE2;
1132 char resname[10];
1133
1134 /* these MUST correspond to the atnms array in do_vsite_aromatics! */
1135 enum { atCG, atND1, atHD1, atCD2, atHD2, atCE1, atHE1, atNE2, atHE2, atNR };
1136 real x[atNR],y[atNR];
1137
1138 /* CG, CE1 and NE2 stay, each gets part of the total mass,
1139 rest gets virtualized */
1140 /* check number of atoms, 3 hydrogens may be missing: */
1141 /* assert( nrfound >= atNR-3 || nrfound <= atNR );
1142 * Don't understand the above logic. Shouldn't it be && rather than || ???
1143 */
1144 if ((nrfound < atNR-3) || (nrfound > atNR))
1145 gmx_incons("Generating vsites for HIS");
1146
1147 /* avoid warnings about uninitialized variables */
1148 b_ND1_HD1=b_NE2_HE2=b_CE1_HE1=b_CD2_HD2=a_NE2_CE1_HE1=
1149 a_NE2_CD2_HD2=a_CE1_ND1_HD1=a_CE1_NE2_HE2=0;
1150
1151 if(ats[atHD1]!=NOTSET) {
1152 if(ats[atHE2]!=NOTSET)
1153 sprintf(resname,"HISH");
1154 else
1155 sprintf(resname,"HISA");
1156 } else {
1157 sprintf(resname,"HISB");
1158 }
1159
1160 /* Get geometry from database */
1161 b_CG_ND1=get_ddb_bond(vsitetop,nvsitetop,resname,"CG","ND1");
1162 b_ND1_CE1=get_ddb_bond(vsitetop,nvsitetop,resname,"ND1","CE1");
1163 b_CE1_NE2=get_ddb_bond(vsitetop,nvsitetop,resname,"CE1","NE2");
1164 b_CG_CD2 =get_ddb_bond(vsitetop,nvsitetop,resname,"CG","CD2");
1165 b_CD2_NE2=get_ddb_bond(vsitetop,nvsitetop,resname,"CD2","NE2");
1166 a_CG_ND1_CE1=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,resname,"CG","ND1","CE1");
1167 a_CG_CD2_NE2=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,resname,"CG","CD2","NE2");
1168 a_ND1_CE1_NE2=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,resname,"ND1","CE1","NE2");
1169 a_CE1_NE2_CD2=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,resname,"CE1","NE2","CD2");
1170
1171 if(ats[atHD1]!=NOTSET) {
1172 b_ND1_HD1=get_ddb_bond(vsitetop,nvsitetop,resname,"ND1","HD1");
1173 a_CE1_ND1_HD1=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,resname,"CE1","ND1","HD1");
1174 }
1175 if(ats[atHE2]!=NOTSET) {
1176 b_NE2_HE2=get_ddb_bond(vsitetop,nvsitetop,resname,"NE2","HE2");
1177 a_CE1_NE2_HE2=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,resname,"CE1","NE2","HE2");
1178 }
1179 if(ats[atHD2]!=NOTSET) {
1180 b_CD2_HD2=get_ddb_bond(vsitetop,nvsitetop,resname,"CD2","HD2");
1181 a_NE2_CD2_HD2=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,resname,"NE2","CD2","HD2");
1182 }
1183 if(ats[atHE1]!=NOTSET) {
1184 b_CE1_HE1=get_ddb_bond(vsitetop,nvsitetop,resname,"CE1","HE1");
1185 a_NE2_CE1_HE1=DEG2RAD*get_ddb_angle(vsitetop,nvsitetop,resname,"NE2","CE1","HE1");
1186 }
1187
1188 /* constraints between CG, CE1 and NE1 */
1189 dCGCE1 = sqrt( cosrule(b_CG_ND1,b_ND1_CE1,a_CG_ND1_CE1) );
1190 dCGNE2 = sqrt( cosrule(b_CG_CD2,b_CD2_NE2,a_CG_CD2_NE2) );
1191
1192 my_add_param(&(plist[F_CONSTRNC]),ats[atCG] ,ats[atCE1],dCGCE1);
1193 my_add_param(&(plist[F_CONSTRNC]),ats[atCG] ,ats[atNE2],dCGNE2);
1194 /* we already have a constraint CE1-NE2, so we don't add it again */
1195
1196 /* calculate the positions in a local frame of reference.
1197 * The x-axis is the line from CG that makes a right angle
1198 * with the bond CE1-NE2, and the y-axis the bond CE1-NE2.
1199 */
1200 /* First calculate the x-axis intersection with y-axis (=yCE1).
1201 * Get cos(angle CG-CE1-NE2) :
1202 */
1203 cosalpha=acosrule(dCGNE2,dCGCE1,b_CE1_NE2);
1204 x[atCE1] = 0;
1205 y[atCE1] = cosalpha*dCGCE1;
1206 x[atNE2] = 0;
1207 y[atNE2] = y[atCE1]-b_CE1_NE2;
1208 sinalpha=sqrt(1-cosalpha*cosalpha);
1209 x[atCG] = - sinalpha*dCGCE1;
1210 y[atCG] = 0;
1211
1212 /* calculate ND1 and CD2 positions from CE1 and NE2 */
1213
1214 x[atND1] = -b_ND1_CE1*sin(a_ND1_CE1_NE2);
1215 y[atND1] = y[atCE1]-b_ND1_CE1*cos(a_ND1_CE1_NE2);
1216
1217 x[atCD2] = -b_CD2_NE2*sin(a_CE1_NE2_CD2);
1218 y[atCD2] = y[atNE2]+b_CD2_NE2*cos(a_CE1_NE2_CD2);
1219
1220 /* And finally the hydrogen positions */
1221 if(ats[atHE1]!=NOTSET) {
1222 x[atHE1] = x[atCE1] + b_CE1_HE1*sin(a_NE2_CE1_HE1);
1223 y[atHE1] = y[atCE1] - b_CE1_HE1*cos(a_NE2_CE1_HE1);
1224 }
1225 /* HD2 - first get (ccw) angle from (positive) y-axis */
1226 if(ats[atHD2]!=NOTSET) {
1227 alpha = a_CE1_NE2_CD2 + M_PI - a_NE2_CD2_HD2;
1228 x[atHD2] = x[atCD2] - b_CD2_HD2*sin(alpha);
1229 y[atHD2] = y[atCD2] + b_CD2_HD2*cos(alpha);
1230 }
1231 if(ats[atHD1]!=NOTSET) {
1232 /* HD1 - first get (cw) angle from (positive) y-axis */
1233 alpha = a_ND1_CE1_NE2 + M_PI - a_CE1_ND1_HD1;
1234 x[atHD1] = x[atND1] - b_ND1_HD1*sin(alpha);
1235 y[atHD1] = y[atND1] - b_ND1_HD1*cos(alpha);
1236 }
1237 if(ats[atHE2]!=NOTSET) {
1238 x[atHE2] = x[atNE2] + b_NE2_HE2*sin(a_CE1_NE2_HE2);
1239 y[atHE2] = y[atNE2] + b_NE2_HE2*cos(a_CE1_NE2_HE2);
1240 }
1241 /* Have all coordinates now */
1242
1243 /* calc center-of-mass; keep atoms CG, CE1, NE2 and
1244 * set the rest to vsite3
1245 */
1246 mtot=xcom=ycom=0;
1247 nvsite=0;
1248 for(i=0; i<atNR; i++)
1249 if (ats[i]!=NOTSET) {
1250 mtot+=at->atom[ats[i]].m;
1251 xcom+=x[i]*at->atom[ats[i]].m;
1252 ycom+=y[i]*at->atom[ats[i]].m;
1253 if (i!=atCG && i!=atCE1 && i!=atNE2) {
1254 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
1255 (*vsite_type)[ats[i]]=F_VSITE3;
1256 nvsite++;
1257 }
1258 }
1259 if (nvsite+3 != nrfound)
1260 gmx_incons("Generating vsites for HIS");
1261
1262 xcom/=mtot;
1263 ycom/=mtot;
1264
1265 /* distribute mass so that com stays the same */
1266 mG = xcom*mtot/x[atCG];
1267 mrest = mtot-mG;
1268 mCE1 = (ycom-y[atNE2])*mrest/(y[atCE1]-y[atNE2]);
1269 mNE2 = mrest-mCE1;
1270
1271 at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = mG;
1272 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB = mCE1;
1273 at->atom[ats[atNE2]].m = at->atom[ats[atNE2]].mB = mNE2;
1274
1275 /* HE1 */
1276 if (ats[atHE1]!=NOTSET) {
1277 calc_vsite3_param(x[atHE1],y[atHE1],x[atCE1],y[atCE1],x[atNE2],y[atNE2],
1278 x[atCG],y[atCG],&a,&b);
1279 add_vsite3_param(&plist[F_VSITE3],
1280 ats[atHE1],ats[atCE1],ats[atNE2],ats[atCG],a,b);
1281 }
1282 /* HE2 */
1283 if (ats[atHE2]!=NOTSET) {
1284 calc_vsite3_param(x[atHE2],y[atHE2],x[atNE2],y[atNE2],x[atCE1],y[atCE1],
1285 x[atCG],y[atCG],&a,&b);
1286 add_vsite3_param(&plist[F_VSITE3],
1287 ats[atHE2],ats[atNE2],ats[atCE1],ats[atCG],a,b);
1288 }
1289
1290 /* ND1 */
1291 calc_vsite3_param(x[atND1],y[atND1],x[atNE2],y[atNE2],x[atCE1],y[atCE1],
1292 x[atCG],y[atCG],&a,&b);
1293 add_vsite3_param(&plist[F_VSITE3],
1294 ats[atND1],ats[atNE2],ats[atCE1],ats[atCG],a,b);
1295
1296 /* CD2 */
1297 calc_vsite3_param(x[atCD2],y[atCD2],x[atCE1],y[atCE1],x[atNE2],y[atNE2],
1298 x[atCG],y[atCG],&a,&b);
1299 add_vsite3_param(&plist[F_VSITE3],
1300 ats[atCD2],ats[atCE1],ats[atNE2],ats[atCG],a,b);
1301
1302 /* HD1 */
1303 if (ats[atHD1]!=NOTSET) {
1304 calc_vsite3_param(x[atHD1],y[atHD1],x[atNE2],y[atNE2],x[atCE1],y[atCE1],
1305 x[atCG],y[atCG],&a,&b);
1306 add_vsite3_param(&plist[F_VSITE3],
1307 ats[atHD1],ats[atNE2],ats[atCE1],ats[atCG],a,b);
1308 }
1309 /* HD2 */
1310 if (ats[atHD2]!=NOTSET) {
1311 calc_vsite3_param(x[atHD2],y[atHD2],x[atCE1],y[atCE1],x[atNE2],y[atNE2],
1312 x[atCG],y[atCG],&a,&b);
1313 add_vsite3_param(&plist[F_VSITE3],
1314 ats[atHD2],ats[atCE1],ats[atNE2],ats[atCG],a,b);
1315 }
1316 return nvsite;
1317 }
1318
1319 static gmx_bool is_vsite(int vsite_type)
1320 {
1321 if (vsite_type == NOTSET)
1322 return FALSE;
1323 switch ( abs(vsite_type) ) {
1324 case F_VSITE3:
1325 case F_VSITE3FD:
1326 case F_VSITE3OUT:
1327 case F_VSITE3FAD:
1328 case F_VSITE4FD:
1329 case F_VSITE4FDN:
1330 return TRUE;
1331 default:
1332 return FALSE;
1333 }
1334 }
1335
1336 static char atomnamesuffix[] = "1234";
1337
1338 void do_vsites(int nrtp, t_restp rtp[], gpp_atomtype_t atype,
1339 t_atoms *at, t_symtab *symtab, rvec *x[],
1340 t_params plist[], int *vsite_type[], int *cgnr[],
1341 real mHmult, gmx_bool bVsiteAromatics,
1342 const char *ffdir)
1343 {
1344 #define MAXATOMSPERRESIDUE 16
1345 int i,j,k,m,i0,ni0,whatres,resind,add_shift,ftype,nvsite,nadd;
1346 int ai,aj,ak,al;
1347 int nrfound=0,needed,nrbonds,nrHatoms,Heavy,nrheavies,tpM,tpHeavy;
1348 int Hatoms[4],heavies[4],bb;
1349 gmx_bool bWARNING,bAddVsiteParam,bFirstWater;
1350 matrix tmpmat;
1351 gmx_bool *bResProcessed;
1352 real mHtot,mtot,fact,fact2;
1353 rvec rpar,rperp,temp;
1354 char name[10],tpname[32],nexttpname[32],*ch;
1355 rvec *newx;
1356 int *o2n,*newvsite_type,*newcgnr,ats[MAXATOMSPERRESIDUE];
1357 t_atom *newatom;
1358 t_params *params;
1359 char ***newatomname;
1360 char *resnm=NULL;
1361 int ndb,f;
1362 char **db;
1363 int nvsiteconf,nvsitetop,cmplength;
1364 gmx_bool isN,planarN,bFound;
1365 gmx_residuetype_t rt;
1366
1367 t_vsiteconf *vsiteconflist;
1368 /* pointer to a list of CH3/NH3/NH2 configuration entries.
1369 * See comments in read_vsite_database. It isnt beautiful,
1370 * but it had to be fixed, and I dont even want to try to
1371 * maintain this part of the code...
1372 */
1373 t_vsitetop *vsitetop;
1374 /* Pointer to a list of geometry (bond/angle) entries for
1375 * residues like PHE, TRP, TYR, HIS, etc., where we need
1376 * to know the geometry to construct vsite aromatics.
1377 * Note that equilibrium geometry isnt necessarily the same
1378 * as the individual bond and angle values given in the
1379 * force field (rings can be strained).
1380 */
1381
1382 /* if bVsiteAromatics=TRUE do_vsites will specifically convert atoms in
1383 PHE, TRP, TYR and HIS to a construction of virtual sites */
1384 enum { resPHE, resTRP, resTYR, resHIS, resNR };
1385 const char *resnms[resNR] = { "PHE", "TRP", "TYR", "HIS" };
1386 /* Amber03 alternative names for termini */
1387 const char *resnmsN[resNR] = { "NPHE", "NTRP", "NTYR", "NHIS" };
1388 const char *resnmsC[resNR] = { "CPHE", "CTRP", "CTYR", "CHIS" };
1389 /* HIS can be known as HISH, HIS1, HISA, HID, HIE, HIP, etc. too */
1390 gmx_bool bPartial[resNR] = { FALSE, FALSE, FALSE, TRUE };
1391 /* the atnms for every residue MUST correspond to the enums in the
1392 gen_vsites_* (one for each residue) routines! */
1393 /* also the atom names in atnms MUST be in the same order as in the .rtp! */
1394 const char *atnms[resNR][MAXATOMSPERRESIDUE+1] = {
1395 { "CG", /* PHE */
1396 "CD1", "HD1", "CD2", "HD2",
1397 "CE1", "HE1", "CE2", "HE2",
1398 "CZ", "HZ", NULL },
1399 { "CB", /* TRP */
1400 "CG",
1401 "CD1", "HD1", "CD2",
1402 "NE1", "HE1", "CE2", "CE3", "HE3",
1403 "CZ2", "HZ2", "CZ3", "HZ3",
1404 "CH2", "HH2", NULL },
1405 { "CG", /* TYR */
1406 "CD1", "HD1", "CD2", "HD2",
1407 "CE1", "HE1", "CE2", "HE2",
1408 "CZ", "OH", "HH", NULL },
1409 { "CG", /* HIS */
1410 "ND1", "HD1", "CD2", "HD2",
1411 "CE1", "HE1", "NE2", "HE2", NULL }
1412 };
1413
1414 if (debug) {
1415 printf("Searching for atoms to make virtual sites ...\n");
1416 fprintf(debug,"# # # VSITES # # #\n");
1417 }
1418
1419 ndb = fflib_search_file_end(ffdir,".vsd",FALSE,&db);
1420 nvsiteconf = 0;
1421 vsiteconflist = NULL;
1422 nvsitetop = 0;
1423 vsitetop = NULL;
1424 for(f=0; f<ndb; f++) {
1425 read_vsite_database(db[f],&vsiteconflist,&nvsiteconf,&vsitetop,&nvsitetop);
1426 sfree(db[f]);
1427 }
1428 sfree(db);
1429
1430 bFirstWater=TRUE;
1431 nvsite=0;
1432 nadd=0;
1433 /* we need a marker for which atoms should *not* be renumbered afterwards */
1434 add_shift = 10*at->nr;
1435 /* make arrays where masses can be inserted into */
1436 snew(newx,at->nr);
1437 snew(newatom,at->nr);
1438 snew(newatomname,at->nr);
1439 snew(newvsite_type,at->nr);
1440 snew(newcgnr,at->nr);
1441 /* make index array to tell where the atoms go to when masses are inserted */
1442 snew(o2n,at->nr);
1443 for(i=0; i<at->nr; i++)
1444 o2n[i]=i;
1445 /* make index to tell which residues were already processed */
1446 snew(bResProcessed,at->nres);
1447
1448 gmx_residuetype_init(&rt);
1449
1450 /* generate vsite constructions */
1451 /* loop over all atoms */
1452 resind = -1;
1453 for(i=0; (i<at->nr); i++) {
1454 if (at->atom[i].resind != resind) {
1455 resind = at->atom[i].resind;
1456 resnm = *(at->resinfo[resind].name);
1457 }
1458 /* first check for aromatics to virtualize */
1459 /* don't waste our effort on DNA, water etc. */
1460 /* Only do the vsite aromatic stuff when we reach the
1461 * CA atom, since there might be an X2/X3 group on the
1462 * N-terminus that must be treated first.
1463 */
1464 if ( bVsiteAromatics &&
1465 !strcmp(*(at->atomname[i]),"CA") &&
1466 !bResProcessed[resind] &&
1467 gmx_residuetype_is_protein(rt,*(at->resinfo[resind].name)) ) {
1468 /* mark this residue */
1469 bResProcessed[resind] = TRUE;
1470 /* find out if this residue needs converting */
1471 whatres=NOTSET;
1472 for(j=0; j<resNR && whatres==NOTSET; j++) {
1473
1474 cmplength = bPartial[j] ? strlen(resnm)-1 : strlen(resnm);
1475
1476 bFound = ((gmx_strncasecmp(resnm,resnms[j], cmplength)==0) ||
1477 (gmx_strncasecmp(resnm,resnmsN[j],cmplength)==0) ||
1478 (gmx_strncasecmp(resnm,resnmsC[j],cmplength)==0));
1479
1480 if ( bFound ) {
1481 whatres=j;
1482 /* get atoms we will be needing for the conversion */
1483 nrfound=0;
1484 for (k=0; atnms[j][k]; k++)
1485 {
1486 ats[k]=NOTSET;
1487 for(m=i; m<at->nr && at->atom[m].resind==resind && ats[k]==NOTSET;m++)
1488 {
1489 if (gmx_strcasecmp(*(at->atomname[m]),atnms[j][k])==0)
1490 {
1491 ats[k]=m;
1492 nrfound++;
1493 }
1494 }
1495 }
1496
1497 /* now k is number of atom names in atnms[j] */
1498 if (j==resHIS)
1499 {
1500 needed = k-3;
1501 }
1502 else
1503 {
1504 needed = k;
1505 }
1506 if (nrfound<needed)
1507 {
1508 gmx_fatal(FARGS,"not enough atoms found (%d, need %d) in "
1509 "residue %s %d while\n "
1510 "generating aromatics virtual site construction",
1511 nrfound,needed,resnm,at->resinfo[resind].nr);
1512 }
1513 /* Advance overall atom counter */
1514 i++;
1515 }
1516 }
1517 /* the enums for every residue MUST correspond to atnms[residue] */
1518 switch (whatres) {
1519 case resPHE:
1520 if (debug) fprintf(stderr,"PHE at %d\n",o2n[ats[0]]+1);
1521 nvsite+=gen_vsites_phe(at, vsite_type, plist, nrfound, ats, vsitetop, nvsitetop);
1522 break;
1523 case resTRP:
1524 if (debug) fprintf(stderr,"TRP at %d\n",o2n[ats[0]]+1);
1525 nvsite+=gen_vsites_trp(atype, &newx, &newatom, &newatomname, &o2n,
1526 &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
1527 at, vsite_type, plist, nrfound, ats, add_shift, vsitetop, nvsitetop);
1528 break;
1529 case resTYR:
1530 if (debug) fprintf(stderr,"TYR at %d\n",o2n[ats[0]]+1);
1531 nvsite+=gen_vsites_tyr(atype, &newx, &newatom, &newatomname, &o2n,
1532 &newvsite_type, &newcgnr, symtab, &nadd, *x, cgnr,
1533 at, vsite_type, plist, nrfound, ats, add_shift, vsitetop, nvsitetop);
1534 break;
1535 case resHIS:
1536 if (debug) fprintf(stderr,"HIS at %d\n",o2n[ats[0]]+1);
1537 nvsite+=gen_vsites_his(at, vsite_type, plist, nrfound, ats, vsitetop, nvsitetop);
1538 break;
1539 case NOTSET:
1540 /* this means this residue won't be processed */
1541 break;
1542 default:
1543 gmx_fatal(FARGS,"DEATH HORROR in do_vsites (%s:%d)",
1544 __FILE__,__LINE__);
1545 } /* switch whatres */
1546 /* skip back to beginning of residue */
1547 while(i>0 && at->atom[i-1].resind == resind)
1548 i--;
1549 } /* if bVsiteAromatics & is protein */
1550
1551 /* now process the rest of the hydrogens */
1552 /* only process hydrogen atoms which are not already set */
1553 if ( ((*vsite_type)[i]==NOTSET) && is_hydrogen(*(at->atomname[i]))) {
1554 /* find heavy atom, count #bonds from it and #H atoms bound to it
1555 and return H atom numbers (Hatoms) and heavy atom numbers (heavies) */
1556 count_bonds(i, &plist[F_BONDS], at->atomname,
1557 &nrbonds, &nrHatoms, Hatoms, &Heavy, &nrheavies, heavies);
1558 /* get Heavy atom type */
1559 tpHeavy=get_atype(Heavy,at,nrtp,rtp,rt);
1560 strcpy(tpname,get_atomtype_name(tpHeavy,atype));
1561
1562 bWARNING=FALSE;
1563 bAddVsiteParam=TRUE;
1564 /* nested if's which check nrHatoms, nrbonds and atomname */
1565 if (nrHatoms == 1) {
1566 switch(nrbonds) {
1567 case 2: /* -O-H */
1568 (*vsite_type)[i]=F_BONDS;
1569 break;
1570 case 3: /* =CH-, -NH- or =NH+- */
1571 (*vsite_type)[i]=F_VSITE3FD;
1572 break;
1573 case 4: /* --CH- (tert) */
1574 /* The old type 4FD had stability issues, so
1575 * all new constructs should use 4FDN
1576 */
1577 (*vsite_type)[i]=F_VSITE4FDN;
1578
1579 /* Check parity of heavy atoms from coordinates */
1580 ai = Heavy;
1581 aj = heavies[0];
1582 ak = heavies[1];
1583 al = heavies[2];
1584 rvec_sub((*x)[aj],(*x)[ai],tmpmat[0]);
1585 rvec_sub((*x)[ak],(*x)[ai],tmpmat[1]);
1586 rvec_sub((*x)[al],(*x)[ai],tmpmat[2]);
1587
1588 if(det(tmpmat)>0)
1589 {
1590 /* swap parity */
1591 heavies[1] = aj;
1592 heavies[0] = ak;
1593 }
1594
1595 break;
1596 default: /* nrbonds != 2, 3 or 4 */
1597 bWARNING=TRUE;
1598 }
1599
1600 } else if ( /*(nrHatoms == 2) && (nrbonds == 2) && REMOVED this test
1601 DvdS 19-01-04 */
1602 (gmx_strncasecmp(*at->atomname[Heavy],"OW",2)==0) ) {
1603 bAddVsiteParam=FALSE; /* this is water: skip these hydrogens */
1604 if (bFirstWater) {
1605 bFirstWater=FALSE;
1606 if (debug)
1607 fprintf(debug,
1608 "Not converting hydrogens in water to virtual sites\n");
1609 }
1610 } else if ( (nrHatoms == 2) && (nrbonds == 4) ) {
1611 /* -CH2- , -NH2+- */
1612 (*vsite_type)[Hatoms[0]] = F_VSITE3OUT;
1613 (*vsite_type)[Hatoms[1]] =-F_VSITE3OUT;
1614 } else {
1615 /* 2 or 3 hydrogen atom, with 3 or 4 bonds in total to the heavy atom.
1616 * If it is a nitrogen, first check if it is planar.
1617 */
1618 isN=planarN=FALSE;
1619 if((nrHatoms == 2) && ((*at->atomname[Heavy])[0]=='N')) {
1620 isN=TRUE;
1621 j=nitrogen_is_planar(vsiteconflist,nvsiteconf,tpname);
1622 if(j<0)
1623 gmx_fatal(FARGS,"No vsite database NH2 entry for type %s\n",tpname);
1624 planarN=(j==1);
1625 }
1626 if ( (nrHatoms == 2) && (nrbonds == 3) && ( !isN || planarN ) ) {
1627 /* =CH2 or, if it is a nitrogen NH2, it is a planar one */
1628 (*vsite_type)[Hatoms[0]] = F_VSITE3FAD;
1629 (*vsite_type)[Hatoms[1]] =-F_VSITE3FAD;
1630 } else if ( ( (nrHatoms == 2) && (nrbonds == 3) &&
1631 ( isN && !planarN ) ) ||
1632 ( (nrHatoms == 3) && (nrbonds == 4) ) ) {
1633 /* CH3, NH3 or non-planar NH2 group */
1634 int Hat_vsite_type[3] = { F_VSITE3, F_VSITE3OUT, F_VSITE3OUT };
1635 gmx_bool Hat_SwapParity[3] = { FALSE, TRUE, FALSE };
1636
1637 if (debug) fprintf(stderr,"-XH3 or nonplanar NH2 group at %d\n",i+1);
1638 bAddVsiteParam=FALSE; /* we'll do this ourselves! */
1639 /* -NH2 (umbrella), -NH3+ or -CH3 */
1640 (*vsite_type)[Heavy] = F_VSITE3;
1641 for (j=0; j<nrHatoms; j++)
1642 (*vsite_type)[Hatoms[j]] = Hat_vsite_type[j];
1643 /* get dummy mass type from first char of heavy atom type (N or C) */
1644
1645 strcpy(nexttpname,get_atomtype_name(get_atype(heavies[0],at,nrtp,rtp,rt),atype));
1646 ch=get_dummymass_name(vsiteconflist,nvsiteconf,tpname,nexttpname);
1647
1648 if (ch == NULL) {
1649 if (ndb > 0) {
1650 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);
1651 } else {
1652 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);
1653 }
1654 } else {
1655 strcpy(name,ch);
1656 }
1657
1658 tpM=vsite_nm2type(name,atype);
1659 /* make space for 2 masses: shift all atoms starting with 'Heavy' */
1660 #define NMASS 2
1661 i0=Heavy;
1662 ni0=i0+nadd;
1663 if (debug)
1664 fprintf(stderr,"Inserting %d dummy masses at %d\n",NMASS,o2n[i0]+1);
1665 nadd+=NMASS;
1666 for(j=i0; j<at->nr; j++)
1667 o2n[j]=j+nadd;
1668
1669 srenew(newx,at->nr+nadd);
1670 srenew(newatom,at->nr+nadd);
1671 srenew(newatomname,at->nr+nadd);
1672 srenew(newvsite_type,at->nr+nadd);
1673 srenew(newcgnr,at->nr+nadd);
1674
1675 for(j=0; j<NMASS; j++)
1676 newatomname[at->nr+nadd-1-j]=NULL;
1677
1678 /* calculate starting position for the masses */
1679 mHtot=0;
1680 /* get atom masses, and set Heavy and Hatoms mass to zero */
1681 for(j=0; j<nrHatoms; j++) {
1682 mHtot += get_amass(Hatoms[j],at,nrtp,rtp,rt);
1683 at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
1684 }
1685 mtot = mHtot + get_amass(Heavy,at,nrtp,rtp,rt);
1686 at->atom[Heavy].m = at->atom[Heavy].mB = 0;
1687 if (mHmult != 1.0)
1688 mHtot *= mHmult;
1689 fact2=mHtot/mtot;
1690 fact=sqrt(fact2);
1691 /* generate vectors parallel and perpendicular to rotational axis:
1692 * rpar = Heavy -> Hcom
1693 * rperp = Hcom -> H1 */
1694 clear_rvec(rpar);
1695 for(j=0; j<nrHatoms; j++)
1696 rvec_inc(rpar,(*x)[Hatoms[j]]);
1697 svmul(1.0/nrHatoms,rpar,rpar); /* rpar = ( H1+H2+H3 ) / 3 */
1698 rvec_dec(rpar,(*x)[Heavy]); /* - Heavy */
1699 rvec_sub((*x)[Hatoms[0]],(*x)[Heavy],rperp);
1700 rvec_dec(rperp,rpar); /* rperp = H1 - Heavy - rpar */
1701 /* calc mass positions */
1702 svmul(fact2,rpar,temp);
1703 for (j=0; (j<NMASS); j++) /* xM = xN + fact2 * rpar +/- fact * rperp */
1704 rvec_add((*x)[Heavy],temp,newx[ni0+j]);
1705 svmul(fact,rperp,temp);
1706 rvec_inc(newx[ni0 ],temp);
1707 rvec_dec(newx[ni0+1],temp);
1708 /* set atom parameters for the masses */
1709 for(j=0; (j<NMASS); j++) {
1710 /* make name: "M??#" or "M?#" (? is atomname, # is number) */
1711 name[0]='M';
1712 for(k=0; (*at->atomname[Heavy])[k] && ( k < NMASS ); k++ )
1713 name[k+1]=(*at->atomname[Heavy])[k];
1714 name[k+1]=atomnamesuffix[j];
1715 name[k+2]='\0';
1716 newatomname[ni0+j] = put_symtab(symtab,name);
1717 newatom[ni0+j].m = newatom[ni0+j].mB = mtot/NMASS;
1718 newatom[ni0+j].q = newatom[ni0+j].qB = 0.0;
1719 newatom[ni0+j].type = newatom[ni0+j].typeB = tpM;
1720 newatom[ni0+j].ptype = eptAtom;
1721 newatom[ni0+j].resind= at->atom[i0].resind;
1722 newvsite_type[ni0+j] = NOTSET;
1723 newcgnr[ni0+j] = (*cgnr)[i0];
1724 }
1725 /* add constraints between dummy masses and to heavies[0] */
1726 /* 'add_shift' says which atoms won't be renumbered afterwards */
1727 my_add_param(&(plist[F_CONSTRNC]),heavies[0], add_shift+ni0, NOTSET);
1728 my_add_param(&(plist[F_CONSTRNC]),heavies[0], add_shift+ni0+1,NOTSET);
1729 my_add_param(&(plist[F_CONSTRNC]),add_shift+ni0,add_shift+ni0+1,NOTSET);
1730
1731 /* generate Heavy, H1, H2 and H3 from M1, M2 and heavies[0] */
1732 /* note that vsite_type cannot be NOTSET, because we just set it */
1733 add_vsite3_atoms (&plist[(*vsite_type)[Heavy]],
1734 Heavy, heavies[0], add_shift+ni0, add_shift+ni0+1,
1735 FALSE);
1736 for(j=0; j<nrHatoms; j++)
1737 add_vsite3_atoms(&plist[(*vsite_type)[Hatoms[j]]],
1738 Hatoms[j], heavies[0], add_shift+ni0, add_shift+ni0+1,
1739 Hat_SwapParity[j]);
1740 #undef NMASS
1741 } else
1742 bWARNING=TRUE;
1743
1744 }
1745 if (bWARNING)
1746 fprintf(stderr,
1747 "Warning: cannot convert atom %d %s (bound to a heavy atom "
1748 "%s with \n"
1749 " %d bonds and %d bound hydrogens atoms) to virtual site\n",
1750 i+1,*(at->atomname[i]),tpname,nrbonds,nrHatoms);
1751 if (bAddVsiteParam) {
1752 /* add vsite parameters to topology,
1753 also get rid of negative vsite_types */
1754 add_vsites(plist, (*vsite_type), Heavy, nrHatoms, Hatoms,
1755 nrheavies, heavies);
1756 /* transfer mass of virtual site to Heavy atom */
1757 for(j=0; j<nrHatoms; j++)
1758 if (is_vsite((*vsite_type)[Hatoms[j]])) {
1759 at->atom[Heavy].m += at->atom[Hatoms[j]].m;
1760 at->atom[Heavy].mB = at->atom[Heavy].m;
1761 at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
1762 }
1763 }
1764 nvsite+=nrHatoms;
1765 if (debug) {
1766 fprintf(debug,"atom %d: ",o2n[i]+1);
1767 print_bonds(debug,o2n,nrHatoms,Hatoms,Heavy,nrheavies,heavies);
1768 }
1769 } /* if vsite NOTSET & is hydrogen */
1770
1771 } /* for i < at->nr */
1772
1773 gmx_residuetype_destroy(rt);
1774
1775 if (debug) {
1776 fprintf(debug,"Before inserting new atoms:\n");
1777 for(i=0; i<at->nr; i++)
1778 fprintf(debug,"%4d %4d %4s %4d %4s %6d %-10s\n",i+1,o2n[i]+1,
1779 at->atomname[i]?*(at->atomname[i]):"(NULL)",
1780 at->resinfo[at->atom[i].resind].nr,
1781 at->resinfo[at->atom[i].resind].name ?
1782 *(at->resinfo[at->atom[i].resind].name):"(NULL)",
1783 (*cgnr)[i],
1784 ((*vsite_type)[i]==NOTSET) ?
1785 "NOTSET" : interaction_function[(*vsite_type)[i]].name);
1786 fprintf(debug,"new atoms to be inserted:\n");
1787 for(i=0; i<at->nr+nadd; i++)
1788 if (newatomname[i])
1789 fprintf(debug,"%4d %4s %4d %6d %-10s\n",i+1,
1790 newatomname[i]?*(newatomname[i]):"(NULL)",
1791 newatom[i].resind,newcgnr[i],
1792 (newvsite_type[i]==NOTSET) ?
1793 "NOTSET" : interaction_function[newvsite_type[i]].name);
1794 }
1795
1796 /* add all original atoms to the new arrays, using o2n index array */
1797 for(i=0; i<at->nr; i++) {
1798 newatomname [o2n[i]] = at->atomname [i];
1799 newatom [o2n[i]] = at->atom [i];
1800 newvsite_type[o2n[i]] = (*vsite_type)[i];
1801 newcgnr [o2n[i]] = (*cgnr) [i];
1802 copy_rvec((*x)[i],newx[o2n[i]]);
1803 }
1804 /* throw away old atoms */
1805 sfree(at->atom);
1806 sfree(at->atomname);
1807 sfree(*vsite_type);
1808 sfree(*cgnr);
1809 sfree(*x);
1810 /* put in the new ones */
1811 at->nr += nadd;
1812 at->atom = newatom;
1813 at->atomname = newatomname;
1814 *vsite_type = newvsite_type;
1815 *cgnr = newcgnr;
1816 *x = newx;
1817 if (at->nr > add_shift)
1818 gmx_fatal(FARGS,"Added impossible amount of dummy masses "
1819 "(%d on a total of %d atoms)\n",nadd,at->nr-nadd);
1820
1821 if (debug) {
1822 fprintf(debug,"After inserting new atoms:\n");
1823 for(i=0; i<at->nr; i++)
1824 fprintf(debug,"%4d %4s %4d %4s %6d %-10s\n",i+1,
1825 at->atomname[i]?*(at->atomname[i]):"(NULL)",
1826 at->resinfo[at->atom[i].resind].nr,
1827 at->resinfo[at->atom[i].resind].name ?
1828 *(at->resinfo[at->atom[i].resind].name):"(NULL)",
1829 (*cgnr)[i],
1830 ((*vsite_type)[i]==NOTSET) ?
1831 "NOTSET" : interaction_function[(*vsite_type)[i]].name);
1832 }
1833
1834 /* now renumber all the interactions because of the added atoms */
1835 for (ftype=0; ftype<F_NRE; ftype++) {
1836 params=&(plist[ftype]);
1837 if (debug)
1838 fprintf(debug,"Renumbering %d %s\n", params->nr,
1839 interaction_function[ftype].longname);
1840 for (i=0; i<params->nr; i++) {
1841 for (j=0; j<NRAL(ftype); j++)
1842 if (params->param[i].a[j]>=add_shift) {
1843 if (debug) fprintf(debug," [%u -> %u]",params->param[i].a[j],
1844 params->param[i].a[j]-add_shift);
1845 params->param[i].a[j]=params->param[i].a[j]-add_shift;
1846 } else {
1847 if (debug) fprintf(debug," [%u -> %d]",params->param[i].a[j],
1848 o2n[params->param[i].a[j]]);
1849 params->param[i].a[j]=o2n[params->param[i].a[j]];
1850 }
1851 if (debug) fprintf(debug,"\n");
1852 }
1853 }
1854 /* now check if atoms in the added constraints are in increasing order */
1855 params=&(plist[F_CONSTRNC]);
1856 for(i=0; i<params->nr; i++)
1857 if ( params->param[i].AI > params->param[i].AJ ) {
1858 j = params->param[i].AJ;
1859 params->param[i].AJ = params->param[i].AI;
1860 params->param[i].AI = j;
1861 }
1862
1863 /* clean up */
1864 sfree(o2n);
1865
1866 /* tell the user what we did */
1867 fprintf(stderr,"Marked %d virtual sites\n",nvsite);
1868 fprintf(stderr,"Added %d dummy masses\n",nadd);
1869 fprintf(stderr,"Added %d new constraints\n",plist[F_CONSTRNC].nr);
1870 }
1871
1872 void do_h_mass(t_params *psb, int vsite_type[], t_atoms *at, real mHmult,
1873 gmx_bool bDeuterate)
1874 {
1875 int i,j,a;
1876
1877 /* loop over all atoms */
1878 for (i=0; i<at->nr; i++)
1879 /* adjust masses if i is hydrogen and not a virtual site */
1880 if ( !is_vsite(vsite_type[i]) && is_hydrogen(*(at->atomname[i])) ) {
1881 /* find bonded heavy atom */
1882 a=NOTSET;
1883 for(j=0; (j<psb->nr) && (a==NOTSET); j++) {
1884 /* if other atom is not a virtual site, it is the one we want */
1885 if ( (psb->param[j].AI==i) &&
1886 !is_vsite(vsite_type[psb->param[j].AJ]) )
1887 a=psb->param[j].AJ;
1888 else if ( (psb->param[j].AJ==i) &&
1889 !is_vsite(vsite_type[psb->param[j].AI]) )
1890 a=psb->param[j].AI;
1891 }
1892 if (a==NOTSET)
1893 gmx_fatal(FARGS,"Unbound hydrogen atom (%d) found while adjusting mass",
1894 i+1);
1895
1896 /* adjust mass of i (hydrogen) with mHmult
1897 and correct mass of a (bonded atom) with same amount */
1898 if (!bDeuterate) {
1899 at->atom[a].m -= (mHmult-1.0)*at->atom[i].m;
1900 at->atom[a].mB-= (mHmult-1.0)*at->atom[i].m;
1901 }
1902 at->atom[i].m *= mHmult;
1903 at->atom[i].mB*= mHmult;
1904 }
1905 }
1906
got merge a long time ago