?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / kernel / gen_dum.c
blob8a537e907c4c69522df2d9f60f51883d4da9988d
1 /*
2 * $Id$
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 2.0
11 *
12 * Copyright (c) 1991-1999
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
15 *
16 * Please refer to:
17 * GROMACS: A message-passing parallel molecular dynamics implementation
18 * H.J.C. Berendsen, D. van der Spoel and R. van Drunen
19 * Comp. Phys. Comm. 91, 43-56 (1995)
20 *
21 * Also check out our WWW page:
22 * http://md.chem.rug.nl/~gmx
23 * or e-mail to:
24 * gromacs@chem.rug.nl
25 *
26 * And Hey:
27 * GRowing Old MAkes el Chrono Sweat
28 */
29 static char *SRCID_gen_dum_c = "$Id$";
30
31 #include <stdio.h>
32 #include <math.h>
33 #include <string.h>
34 #include "assert.h"
35 #include "gen_dum.h"
36 #include "smalloc.h"
37 #include "resall.h"
38 #include "add_par.h"
39 #include "vec.h"
40 #include "toputil.h"
41 #include "physics.h"
42 #include "index.h"
43 #include "names.h"
44
45 static void count_bonds(int atom, t_params *psb, char ***atomname,
46 int *nrbonds, int *nrHatoms, int Hatoms[], int *Heavy,
47 int *nrheavies, int heavies[])
48 {
49 int i,heavy,other,nrb,nrH,nrhv;
50
51 /* find heavy atom bound to this hydrogen */
52 heavy=NOTSET;
53 for(i=0; (i<psb->nr) && (heavy==NOTSET); i++)
54 if (psb->param[i].AI==atom)
55 heavy=psb->param[i].AJ;
56 else if (psb->param[i].AJ==atom)
57 heavy=psb->param[i].AI;
58 if (heavy==NOTSET)
59 fatal_error(0,"unbound hydrogen atom %d",atom+1);
60 /* find all atoms bound to heavy atom */
61 other=NOTSET;
62 nrb=0;
63 nrH=0;
64 nrhv=0;
65 for(i=0; i<psb->nr; i++) {
66 if (psb->param[i].AI==heavy)
67 other=psb->param[i].AJ;
68 else if (psb->param[i].AJ==heavy)
69 other=psb->param[i].AI;
70 if (other!=NOTSET) {
71 nrb++;
72 if (is_hydrogen(*(atomname[other]))) {
73 Hatoms[nrH]=other;
74 nrH++;
75 } else {
76 heavies[nrhv]=other;
77 nrhv++;
78 }
79 other=NOTSET;
80 }
81 }
82 *Heavy =heavy;
83 *nrbonds =nrb;
84 *nrHatoms=nrH;
85 *nrheavies=nrhv;
86 }
87
88 static void print_bonds(FILE *fp, int o2n[],
89 int nrHatoms, int Hatoms[], int Heavy,
90 int nrheavies, int heavies[])
91 {
92 int i;
93
94 fprintf(fp,"Found: %d Hatoms: ",nrHatoms);
95 for(i=0; i<nrHatoms; i++)
96 fprintf(fp," %d",o2n[Hatoms[i]]+1);
97 fprintf(fp,"; %d Heavy atoms: %d",nrheavies+1,o2n[Heavy]+1);
98 for(i=0; i<nrheavies; i++)
99 fprintf(fp," %d",o2n[heavies[i]]+1);
100 fprintf(fp,"\n");
101 }
102
103 static int get_atype(int atom, t_atoms *at, int nrtp, t_restp rtp[])
104 {
105 int type;
106 bool bNterm;
107 int j;
108 t_restp *rtpp;
109
110 if (at->atom[atom].m)
111 type=at->atom[atom].type;
112 else {
113 /* get type from rtp */
114 rtpp = search_rtp(*(at->resname[at->atom[atom].resnr]),nrtp,rtp);
115 bNterm = is_protein(*(at->resname[at->atom[atom].resnr])) &&
116 (at->atom[atom].resnr == 0);
117 j=search_jtype(rtpp,*(at->atomname[atom]),bNterm);
118 type=rtpp->atom[j].type;
119 }
120 return type;
121 }
122
123 static int nm2type(char *name, t_atomtype *atype)
124 {
125 int tp,j;
126
127 tp=NOTSET;
128 for(j=0; (j < atype->nr) && (tp==NOTSET); j++)
129 if (strcasecmp(name,*(atype->atomname[j])) == 0)
130 tp=j;
131 if (tp==NOTSET)
132 fatal_error(0,"Dummy mass type (%s) not found in atom type database",name);
133 return tp;
134 }
135
136 static real get_amass(int atom, t_atoms *at, int nrtp, t_restp rtp[])
137 {
138 real mass;
139 bool bNterm;
140 int j;
141 t_restp *rtpp;
142
143 if (at->atom[atom].m)
144 mass=at->atom[atom].m;
145 else {
146 /* get mass from rtp */
147 rtpp = search_rtp(*(at->resname[at->atom[atom].resnr]),nrtp,rtp);
148 bNterm = is_protein(*(at->resname[at->atom[atom].resnr])) &&
149 (at->atom[atom].resnr == 0);
150 j=search_jtype(rtpp,*(at->atomname[atom]),bNterm);
151 mass=rtpp->atom[j].m;
152 }
153 return mass;
154 }
155
156 static void my_add_param(t_params *plist, int ai, int aj, real b)
157 {
158 static real c[MAXFORCEPARAM] =
159 { NOTSET, NOTSET, NOTSET, NOTSET, NOTSET, NOTSET };
160
161 c[0]=b;
162 add_param(plist,ai,aj,c,NULL);
163 }
164
165 static void add_dum_atoms(t_params plist[], int dummy_type[],
166 int Heavy, int nrHatoms, int Hatoms[],
167 int nrheavies, int heavies[])
168 {
169 int i,j,ftype,other,moreheavy,bb;
170 bool bSwapParity;
171
172 for(i=0; i<nrHatoms; i++) {
173 ftype=dummy_type[Hatoms[i]];
174 bSwapParity = (ftype<0);
175 dummy_type[Hatoms[i]] = ftype = abs(ftype);
176 if (ftype == F_BONDS) {
177 if ( (nrheavies != 1) && (nrHatoms != 1) )
178 fatal_error(0,"cannot make constraint in add_dum_atoms for %d heavy "
179 "atoms and %d hydrogen atoms",nrheavies,nrHatoms);
180 my_add_param(&(plist[F_SHAKENC]),Hatoms[i],heavies[0],NOTSET);
181 } else {
182 switch (ftype) {
183 case F_DUMMY3:
184 case F_DUMMY3FD:
185 case F_DUMMY3OUT:
186 if (nrheavies < 2)
187 fatal_error(0,"Not enough heavy atoms (%d) for %s (min 3)",
188 nrheavies+1,
189 interaction_function[dummy_type[Hatoms[i]]].name);
190 add_dum3_atoms(&plist[ftype],Hatoms[i],Heavy,heavies[0],heavies[1],
191 bSwapParity);
192 break;
193 case F_DUMMY3FAD: {
194 if (nrheavies > 1)
195 moreheavy=heavies[1];
196 else {
197 /* find more heavy atoms */
198 other=moreheavy=NOTSET;
199 for(j=0; (j<plist[F_BONDS].nr) && (moreheavy==NOTSET); j++) {
200 if (plist[F_BONDS].param[j].AI==heavies[0])
201 other=plist[F_BONDS].param[j].AJ;
202 else if (plist[F_BONDS].param[j].AJ==heavies[0])
203 other=plist[F_BONDS].param[j].AI;
204 if ( (other != NOTSET) && (other != Heavy) )
205 moreheavy=other;
206 }
207 if (moreheavy==NOTSET)
208 fatal_error(0,"Unbound molecule part %d-%d",Heavy+1,Hatoms[0]+1);
209 }
210 add_dum3_atoms(&plist[ftype],Hatoms[i],Heavy,heavies[0],moreheavy,
211 bSwapParity);
212 break;
213 }
214 case F_DUMMY4FD: {
215 if (nrheavies < 3)
216 fatal_error(0,"Not enough heavy atoms (%d) for %s (min 4)",
217 nrheavies+1,
218 interaction_function[dummy_type[Hatoms[i]]].name);
219 add_dum4_atoms(&plist[ftype],
220 Hatoms[0], Heavy, heavies[0], heavies[1], heavies[2]);
221 break;
222 }
223 default:
224 fatal_error(0,"can't use add_dum_atoms for interaction function %s",
225 interaction_function[dummy_type[Hatoms[i]]].name);
226 } /* switch ftype */
227 } /* else */
228 } /* for i */
229 }
230
231 /* some things used by gen_dums_*
232 these are the default GROMACS bondlengths and angles for aromatic
233 ring systems, which are also identical to the GROMOS numbers.
234 Probably these should be in a nice file for easy editing. */
235 #define bR6 0.139
236 #define bR5 0.133
237 #define bCC 0.139
238 #define bCO 0.136
239 #define bCH 0.108
240 #define bNH 0.100
241 #define bOH 0.100
242 /* it does not make sense to change these angles: */
243 #define aR6 (DEG2RAD*120)
244 #define aR6H (M_PI-0.5*aR6)
245 #define aR5 (DEG2RAD*108)
246 #define aR5H (M_PI-0.5*aR5)
247 #define aCOH (DEG2RAD*109.5)
248 /* cosine rule: a^2 = b^2 + c^2 - 2 b c cos(alpha) */
249 /* get a^2 when a, b and alpha are given: */
250 #define cosrule(b,c,alpha) ( sqr(b) + sqr(c) - 2*b*c*cos(alpha) )
251 /* get cos(alpha) when a, b and c are given: */
252 #define acosrule(a,b,c) ( (sqr(b)+sqr(c)-sqr(a))/(2*b*c) )
253
254 static int gen_dums_6ring(t_atoms *at, int *dummy_type[], t_params plist[],
255 int nrfound, int *ats, bool bDoZ)
256 {
257 /* these MUST correspond to the atnms array in do_dum_aromatics! */
258 enum { atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
259 atCZ, atHZ, atNR };
260
261 int i,ndum;
262 real a,b,dCGCE,tmp1,tmp2,mtot;
263 /* CG, CE1 and CE2 stay and each get 1/3 of the total mass,
264 rest gets dummified */
265
266 if (bDoZ) {
267 assert(atNR == nrfound);
268 }
269
270 /* constraints between CG, CE1 and CE2: */
271 dCGCE = sqrt( cosrule(bR6,bR6,aR6) );
272 my_add_param(&(plist[F_SHAKENC]),ats[atCG] ,ats[atCE1],dCGCE);
273 my_add_param(&(plist[F_SHAKENC]),ats[atCG] ,ats[atCE2],dCGCE);
274 my_add_param(&(plist[F_SHAKENC]),ats[atCE1],ats[atCE2],dCGCE);
275
276 /* rest will be dummy3 */
277 mtot=0;
278 ndum=0;
279 for(i=0; i<atNR; i++) {
280 mtot+=at->atom[ats[i]].m;
281 if ( i!=atCG && i!=atCE1 && i!=atCE2 && (bDoZ || (i!=atHZ && i!=atCZ) ) ) {
282 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
283 (*dummy_type)[ats[i]]=F_DUMMY3;
284 ndum++;
285 }
286 }
287 at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB =
288 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB =
289 at->atom[ats[atCE2]].m = at->atom[ats[atCE2]].mB = mtot / 3;
290
291 /* dummy3 construction: r_d = r_i + a r_ij + b r_ik */
292 tmp1 = dCGCE*sin(DEG2RAD*60);
293 tmp2 = bR6*cos(0.5*aR6) + tmp1;
294 tmp1 *= 2;
295 a = b = - bCH / tmp1;
296 /* HE1 and HE2: */
297 add_dum3_param(&plist[F_DUMMY3],
298 ats[atHE1],ats[atCE1],ats[atCE2],ats[atCG], a,b);
299 add_dum3_param(&plist[F_DUMMY3],
300 ats[atHE2],ats[atCE2],ats[atCE1],ats[atCG], a,b);
301 /* CD1, CD2 and CZ: */
302 a = b = tmp2 / tmp1;
303 add_dum3_param(&plist[F_DUMMY3],
304 ats[atCD1],ats[atCE2],ats[atCE1],ats[atCG], a,b);
305 add_dum3_param(&plist[F_DUMMY3],
306 ats[atCD2],ats[atCE1],ats[atCE2],ats[atCG], a,b);
307 if (bDoZ)
308 add_dum3_param(&plist[F_DUMMY3],
309 ats[atCZ], ats[atCG], ats[atCE1],ats[atCE2],a,b);
310 /* HD1, HD2 and HZ: */
311 a = b = ( bCH + tmp2 ) / tmp1;
312 add_dum3_param(&plist[F_DUMMY3],
313 ats[atHD1],ats[atCE2],ats[atCE1],ats[atCG], a,b);
314 add_dum3_param(&plist[F_DUMMY3],
315 ats[atHD2],ats[atCE1],ats[atCE2],ats[atCG], a,b);
316 if (bDoZ)
317 add_dum3_param(&plist[F_DUMMY3],
318 ats[atHZ], ats[atCG], ats[atCE1],ats[atCE2],a,b);
319
320 return ndum;
321 }
322
323 static int gen_dums_phe(t_atoms *at, int *dummy_type[], t_params plist[],
324 int nrfound, int *ats)
325 {
326 return gen_dums_6ring(at, dummy_type, plist, nrfound, ats, TRUE);
327 }
328
329 static int gen_dums_trp(t_atomtype *atype, rvec *newx[],
330 t_atom *newatom[], char ***newatomname[],
331 int *o2n[], int *newdummy_type[], int *newcgnr[],
332 t_symtab *symtab, int *nadd, rvec x[], int *cgnr[],
333 t_atoms *at, int *dummy_type[], t_params plist[],
334 int nrfound, int *ats, int add_shift)
335 {
336 #define NMASS 2
337 /* these MUST correspond to the atnms array in do_dum_aromatics! */
338 enum {
339 atCB, atCG, atCD1, atHD1, atCD2, atNE1, atHE1, atCE2, atCE3, atHE3,
340 atCZ2, atHZ2, atCZ3, atHZ3, atCH2, atHH2, atNR };
341 /* weights for determining the COM's of both rings (M1 and M2): */
342 real mw[NMASS][atNR] = {
343 { 0, 1, 1, 1, 0.5, 1, 1, 0.5, 0, 0,
344 0, 0, 0, 0, 0, 0 },
345 { 0, 0, 0, 0, 0.5, 0, 0, 0.5, 1, 1,
346 1, 1, 1, 1, 1, 1 }
347 };
348 /* flag which sets atoms to positive or negative y-axis: */
349 int py[atNR] = {
350 1, 1, 1, 1, 1, -1, -1, -1, 1, 1,
351 -1, -1, 1, 1, -1, -1
352 };
353 real xi[atNR],yi[atNR],xM[NMASS];
354 int atM[NMASS],tpM,i,i0,j,ndum;
355 real mwtot,mtot,mM[NMASS],xcom,dCBM1,dCBM2,dM1M2;
356 real a,b,c[MAXFORCEPARAM];
357 rvec rx,rcom,temp;
358 char name[10];
359
360 /* keep CB and generate two dummy masses (one in each ring) on the
361 symmetry axis (line through CD1 and middle of CZ3-CH bond) such that
362 center of mass (COM) remains the same and the dummy mass for the
363 five-ring is perpendicularly placed to the CG atom. This results
364 in a somewhat larger moment of inertia, but gives two perpendicular
365 construction vectors for the dummy atoms which makes calculating
366 the dummy constructions easier. */
367
368 assert(atNR == nrfound);
369
370 /* get dummy mass type */
371 tpM=nm2type("MW",atype);
372 /* make space for 2 masses: shift all atoms starting with CB */
373 i0=ats[atCB];
374 for(j=0; j<NMASS; j++)
375 atM[j]=i0+*nadd+j;
376 if (debug)
377 fprintf(stderr,"Inserting %d dummy masses at %d\n",NMASS,(*o2n)[i0]+1);
378 *nadd+=NMASS;
379 for(j=i0; j<at->nr; j++)
380 (*o2n)[j]=j+*nadd;
381 srenew(*newx,at->nr+*nadd);
382 srenew(*newatom,at->nr+*nadd);
383 srenew(*newatomname,at->nr+*nadd);
384 srenew(*newdummy_type,at->nr+*nadd);
385 srenew(*newcgnr,at->nr+*nadd);
386 for(j=0; j<NMASS; j++)
387 newatomname[at->nr+*nadd-1-j]=NULL;
388
389 /* calculate the relative positions of all atoms (along symmetry axis): */
390 /* take x=0 at CD2/CE2: */
391 xi[atCD2] = xi[atCE2] = 0;
392 /* first 5-ring: */
393 xi[atCG] = xi[atNE1] = xi[atCD2] - sin(aR5)*bR5;
394 xi[atHE1] = xi[atNE1] - sin(aR5H-aR5)*bNH;
395 xi[atCB] = xi[atCG] - sin(aR5H-aR5)*bCC;
396 xi[atCD1] = xi[atNE1] - sin(2*aR5-M_PI)*bR5;
397 xi[atHD1] = xi[atCD1] - bCH;
398 /* then 6-ring: */
399 xi[atCE3] = xi[atHE3] =
400 xi[atCZ2] = xi[atHZ2] = xi[atCD2] + sin(aR6)*bR6;
401 xi[atCZ3] = xi[atCH2] = xi[atCE3] + sin(2*aR6-M_PI)*bR6;
402 xi[atHZ3] = xi[atHH2] = xi[atCZ3] + sin(aR6H+M_PI-2*aR6)*bCH;
403 /* also position perp. to symmetry axis: */
404 yi[atCD1] = yi[atHD1] = 0;
405 yi[atCG] = yi[atNE1] = sin(0.5*aR5)*bR5;
406 yi[atCB] = yi[atCG] + cos(aR5H-aR5)*bCC;
407 yi[atHE1] = yi[atNE1] + cos(aR5H-aR5)*bNH;
408 yi[atCD2] = yi[atCE2] =
409 yi[atCZ3] = yi[atCH2] = 0.5*bR6;
410 yi[atCE3] = yi[atCZ2] = yi[atCD2] - cos(aR6)*bR6;
411 yi[atHE3] = yi[atHZ2] = yi[atCE3] + cos(aR6H-aR6)*bCH;
412 yi[atHZ3] = yi[atHH2] = yi[atCZ3] - cos(aR6)*bCH;
413 /* now introduce positive or negative y: */
414 for(i=0; i<atNR; i++)
415 yi[i]*=py[i];
416
417 /* first get COM: */
418 xcom=0;
419 for(j=0; j<NMASS; j++)
420 mM[j]=0;
421 for(i=0; i<atNR; i++)
422 if (i!=atCB) {
423 for(j=0; j<NMASS; j++)
424 mM[j] += mw[j][i] * at->atom[ats[i]].m;
425 xcom += at->atom[ats[i]].m * xi[i];
426 }
427 mtot=0;
428 for(j=0; j<NMASS; j++)
429 mtot += mM[j];
430 xcom /= mtot;
431 /* redefine coordinates to get zero at COM: */
432 for(i=0; i<atNR; i++)
433 xi[i] -= xcom;
434 /* now we set M1 at the same 'x' as CB: */
435 xM[0] = xi[atCB];
436 /* then set M2 so that com is conserved (mM1 * xM1 + mM2 * xM2 = 0): */
437 xM[1] = - mM[0] * xM[0] / mM[1];
438
439 /* make a unitvector that defines our symmetry axis: */
440 rvec_sub(x[ats[atCZ3]],x[ats[atCD2]],rx);
441 unitv(rx,rx); /* rx = ( CZ3 - CD2 ) / | CZ3 - CD2 | */
442 /* make vector that points to origin (= COM): */
443 rvec_add(x[ats[atCE2]],x[ats[atCD2]],rcom);
444 svmul(HALF,rcom,rcom);
445 svmul(xcom,rx,temp);
446 rvec_inc(rcom,temp); /* rcom = 0.5 * ( CE2 + CD2 ) + xcom * rx */
447
448 /* calc initial position for dummy masses: */
449 for(j=0; j<NMASS; j++) {
450 svmul(xM[j],rx,temp);
451 rvec_add(rcom,temp,(*newx)[atM[j]]); /* rM = rcom + xM * rx */
452 }
453
454 /* set parameters for the masses */
455 for(j=0; j<NMASS; j++) {
456 sprintf(name,"MW%d",j+1);
457 (*newatomname) [atM[j]] = put_symtab(symtab,name);
458 (*newatom) [atM[j]].m = (*newatom)[atM[j]].mB = mM[j];
459 (*newatom) [atM[j]].q = (*newatom)[atM[j]].qB = 0.0;
460 (*newatom) [atM[j]].type = (*newatom)[atM[j]].typeB = tpM;
461 (*newatom) [atM[j]].ptype = eptAtom;
462 (*newatom) [atM[j]].resnr = at->atom[i0].resnr;
463 (*newatom) [atM[j]].chain = at->atom[i0].chain;
464 (*newdummy_type)[atM[j]] = NOTSET;
465 (*newcgnr) [atM[j]] = (*cgnr)[i0];
466 }
467 /* renumber cgnr: */
468 for(i=i0; i<at->nr; i++)
469 (*cgnr)[i]++;
470
471 /* constraints between CB, M1 and M2 */
472 /* 'add_shift' says which atoms won't be renumbered afterwards */
473 dCBM1 = yi[atCB];
474 dM1M2 = xM[1]-xM[0];
475 dCBM2 = sqrt( sqr(dCBM1) + sqr(dM1M2) );
476 my_add_param(&(plist[F_SHAKENC]), ats[atCB], add_shift+atM[0], dCBM1);
477 my_add_param(&(plist[F_SHAKENC]), ats[atCB], add_shift+atM[1], dCBM2);
478 my_add_param(&(plist[F_SHAKENC]), add_shift+atM[0], add_shift+atM[1], dM1M2);
479
480 /* rest will be dummy3 */
481 ndum=0;
482 for(i=0; i<atNR; i++)
483 if (i!=atCB) {
484 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
485 (*dummy_type)[ats[i]] = F_DUMMY3;
486 ndum++;
487 }
488
489 /* now define all dummies from M1, M2, CB, ie:
490 r_d = r_M1 + a r_M!M2 + b r_M1_CB */
491 for(i=0; i<atNR; i++)
492 if ( (*dummy_type)[ats[i]] == F_DUMMY3 )
493 add_dum3_param(&plist[F_DUMMY3],
494 ats[i],add_shift+atM[0],add_shift+atM[1],ats[atCB],
495 (xi[i]-xM[0])/dM1M2,yi[i]/dCBM1);
496
497 return ndum;
498 #undef NMASS
499 }
500
501 static int gen_dums_tyr(t_atoms *at, int *dummy_type[], t_params plist[],
502 int nrfound, int *ats)
503 {
504 int ndum;
505 real dCGCE,dCEOH,dCGHH,tmp1,a,b;
506
507 /* these MUST correspond to the atnms array in do_dum_aromatics! */
508 enum { atCG, atCD1, atHD1, atCD2, atHD2, atCE1, atHE1, atCE2, atHE2,
509 atCZ, atOH, atHH, atNR };
510 /* CG, CE1, CE2 (as in general 6-ring) and OH and HH stay,
511 rest gets dummified.
512 Now we have two linked triangles with one improper keeping them flat */
513 assert(atNR == nrfound);
514
515 /* first do 6 ring as default,
516 except CZ (we'll do that different) and HZ (we don't have that): */
517 ndum = gen_dums_6ring(at, dummy_type, plist, nrfound, ats, FALSE);
518
519 /* then construct CZ from the 2nd triangle */
520 /* dummy3 construction: r_d = r_i + a r_ij + b r_ik */
521 a = b = 0.5 * bCO / ( bCO - bR6*cos(aR6) );
522 add_dum3_param(&plist[F_DUMMY3],
523 ats[atCZ],ats[atOH],ats[atCE1],ats[atCE2],a,b);
524 at->atom[ats[atCZ]].m = at->atom[ats[atCZ]].mB = 0;
525
526 /* constraints between CE1, CE2 and OH */
527 dCGCE = sqrt( cosrule(bR6,bR6,aR6) );
528 dCEOH = sqrt( cosrule(bR6,bCO,aR6) );
529 my_add_param(&(plist[F_SHAKENC]),ats[atCE1],ats[atOH],dCEOH);
530 my_add_param(&(plist[F_SHAKENC]),ats[atCE2],ats[atOH],dCEOH);
531 /* assume we also want the COH angle constrained: */
532 tmp1 = bR6*cos(0.5*aR6) + dCGCE*sin(DEG2RAD*60) + bCO;
533 dCGHH = sqrt( cosrule(tmp1,bOH,aCOH) );
534 (*dummy_type)[ats[atHH]]=F_SHAKE;
535 my_add_param(&(plist[F_SHAKENC]),ats[atCG],ats[atHH],dCGHH);
536
537 return ndum;
538 }
539
540 static int gen_dums_his(t_atoms *at, int *dummy_type[], t_params plist[],
541 int nrfound, int *ats)
542 {
543 int ndum,i;
544 real a,b,alpha,dGE,dCENE,mtot;
545 real xCG, yE, mE, mG;
546 real xG, xD, yD, xE, xHD, yHD, xHE, yHE;
547
548 /* these MUST correspond to the atnms array in do_dum_aromatics! */
549 enum { atCG, atND1, atHD1, atCD2, atHD2, atCE1, atHE1, atNE2, atHE2, atNR };
550 /* CG, CE1 and NE2 stay, each gets part of the total mass,
551 rest gets dummified */
552 /* check number of atoms, 3 hydrogens may be missing: */
553 assert( nrfound >= atNR-3 || nrfound <= atNR );
554
555 /* constraints between CG, CE1 and NE1 */
556 dGE = sqrt( cosrule(bR5,bR5,aR5) );
557 dCENE = bR5;
558 my_add_param(&(plist[F_SHAKENC]),ats[atCG] ,ats[atCE1],dGE);
559 my_add_param(&(plist[F_SHAKENC]),ats[atCG] ,ats[atNE2],dGE);
560 /* we already have a constraint CE-NE, so we don't add it again */
561
562 /* rest will be dummy3 */
563 mtot=0;
564 ndum=0;
565 for(i=0; i<atNR; i++)
566 if (ats[i]!=NOTSET) {
567 mtot+=at->atom[ats[i]].m;
568 if (i!=atCG && i!=atCE1 && i!=atNE2) {
569 at->atom[ats[i]].m = at->atom[ats[i]].mB = 0;
570 (*dummy_type)[ats[i]]=F_DUMMY3;
571 ndum++;
572 }
573 }
574 assert(ndum+3 == nrfound);
575
576 /* distribute mass so that com stays the same */
577 xG = bR5 * sin(0.5*aR5) / sin(aR5);
578 xE = xG * sin(0.5*aR5);
579 mE = mtot * xG / ( 2 * ( xG + xE ) );
580 mG = mtot - 2 * mE;
581 at->atom[ats[atCG]].m = at->atom[ats[atCG]].mB = mG;
582 at->atom[ats[atCE1]].m = at->atom[ats[atCE1]].mB =
583 at->atom[ats[atNE2]].m = at->atom[ats[atNE2]].mB = mE;
584
585 /* these we need for all the atoms: */
586 alpha = acos( acosrule(dGE,dGE,dCENE) ); /* angle CG-NE2-CE1 */
587 /* we define 'x' and 'y' perp to the three construction vectors: */
588 xCG = sin(alpha)*dGE; /* perp to r_CE1-NE2 */
589 yE = sin(alpha)*dCENE; /* perp to r_CG-CE1 or r_CG_NE2 */
590
591 /* HE2 */
592 yHE = sin(M_2PI-alpha-aR5H)*bNH;
593 xHE = sin(aR5H)*bNH;
594 a = - yHE / yE;
595 b = - xHE / xCG;
596 if (ats[atHE1]!=NOTSET)
597 add_dum3_param(&plist[F_DUMMY3],
598 ats[atHE1],ats[atCE1],ats[atNE2],ats[atCG],a,b);
599 if (ats[atHE2]!=NOTSET)
600 add_dum3_param(&plist[F_DUMMY3],
601 ats[atHE2],ats[atNE2],ats[atCE1],ats[atCG],a,b);
602
603 /* ND1, CD2 */
604 xD = sin(aR5)*bR5;
605 yD = yE + sin(alpha+aR5-M_PI)*bR5;
606 a = yD / yE;
607 b = xD / xCG;
608 add_dum3_param(&plist[F_DUMMY3],
609 ats[atND1],ats[atNE2],ats[atCE1],ats[atCG],a,b);
610 add_dum3_param(&plist[F_DUMMY3],
611 ats[atCD2],ats[atCE1],ats[atNE2],ats[atCG],a,b);
612
613 /* HD1 */
614 xHD = xD + sin(aR5H-aR5)*bNH;
615 yHD = yD + sin(alpha+aR5-aR5H)*bNH;
616 a = yHD / yE;
617 b = xHD / xCG;
618 if (ats[atHD1]!=NOTSET)
619 add_dum3_param(&plist[F_DUMMY3],
620 ats[atHD1],ats[atNE2],ats[atCE1],ats[atCG],a,b);
621 if (ats[atHD2]!=NOTSET)
622 add_dum3_param(&plist[F_DUMMY3],
623 ats[atHD2],ats[atCE1],ats[atNE2],ats[atCG],a,b);
624
625 return ndum;
626 }
627
628 static bool is_dum(int dummy_type)
629 {
630 if (dummy_type == NOTSET)
631 return FALSE;
632 switch ( abs(dummy_type) ) {
633 case F_DUMMY3:
634 case F_DUMMY3FD:
635 case F_DUMMY3OUT:
636 case F_DUMMY3FAD:
637 case F_DUMMY4FD:
638 return TRUE;
639 default:
640 return FALSE;
641 }
642 }
643
644 /* this seems overdone, but it is FOOLPROOF!!! */
645 static char atomnamesuffix[] = "1234";
646
647 void do_dummies(int nrtp, t_restp rtp[], t_atomtype *atype,
648 t_atoms *at, t_symtab *symtab, rvec *x[],
649 t_params plist[], int *dummy_type[], int *cgnr[],
650 real mHmult, bool bDummyAromatics)
651 {
652 #define MAXATOMSPERRESIDUE 16
653 int i,j,k,i0,ni0,whatres,resnr,add_shift,ftype,ndum,nadd;
654 int nrfound=0,needed,nrbonds,nrHatoms,Heavy,nrheavies,tpM,tpHeavy;
655 int Hatoms[4],heavies[4],bb;
656 bool bWARNING,bAddDumParam,bFirstWater;
657 bool *bResProcessed;
658 real mHtot,mtot,fact,fact2;
659 rvec rpar,rperp,temp;
660 char name[10],tpname[10];
661 rvec *newx;
662 int *o2n,*newdummy_type,*newcgnr,ats[MAXATOMSPERRESIDUE];
663 t_atom *newatom;
664 t_params *params;
665 char ***newatomname;
666 char *resnm=NULL;
667
668 /* if bDummyAromatics=TRUE do_dummies will specifically convert atoms in
669 PHE, TRP, TYR and HIS to a construction of dummy atoms */
670 enum { resPHE, resTRP, resTYR, resHIS, resNR };
671 char *resnms[resNR] = { "PHE", "TRP", "TYR", "HIS" };
672 bool bPartial[resNR] = { FALSE, FALSE, FALSE, TRUE };
673 /* the atnms for every residue MUST correspond to the enums in the
674 gen_dums_* (one for each residue) routines! */
675 /* also the atom names in atnms MUST be in the same order as in the .rtp! */
676 char *atnms[resNR][MAXATOMSPERRESIDUE+1] = {
677 { "CG", /* PHE */
678 "CD1", "HD1", "CD2", "HD2",
679 "CE1", "HE1", "CE2", "HE2",
680 "CZ", "HZ", NULL },
681 { "CB", /* TRP */
682 "CG",
683 "CD1", "HD1", "CD2",
684 "NE1", "HE1", "CE2", "CE3", "HE3",
685 "CZ2", "HZ2", "CZ3", "HZ3",
686 "CH2", "HH2", NULL },
687 { "CG", /* TYR */
688 "CD1", "HD1", "CD2", "HD2",
689 "CE1", "HE1", "CE2", "HE2",
690 "CZ", "OH", "HH", NULL },
691 { "CG", /* HIS */
692 "ND1", "HD1", "CD2", "HD2",
693 "CE1", "HE1", "NE2", "HE2", NULL }
694 };
695
696 if (debug) {
697 printf("Searching for atoms to make dumies...\n");
698 fprintf(debug,"# # # DUMMIES # # #\n");
699 }
700
701 bFirstWater=TRUE;
702 ndum=0;
703 nadd=0;
704 /* we need a marker for which atoms should *not* be renumbered afterwards */
705 add_shift = 10*at->nr;
706 /* make arrays where masses can be inserted into */
707 snew(newx,at->nr);
708 snew(newatom,at->nr);
709 snew(newatomname,at->nr);
710 snew(newdummy_type,at->nr);
711 snew(newcgnr,at->nr);
712 /* make index array to tell where the atoms go to when masses are inse'd */
713 snew(o2n,at->nr);
714 for(i=0; i<at->nr; i++)
715 o2n[i]=i;
716 /* make index to tell which residues were already processed */
717 snew(bResProcessed,at->nres);
718
719 /* generate dummy constructions */
720 /* loop over all atoms */
721 resnr=NOTSET;
722 for(i=0; (i<at->nr); i++) {
723 if (at->atom[i].resnr != resnr) {
724 resnr=at->atom[i].resnr;
725 resnm=*(at->resname[resnr]);
726 }
727 /* first check for aromatics to dummify */
728 /* don't waste our effort on DNA, water etc. */
729 if ( bDummyAromatics &&
730 !bResProcessed[resnr] && is_protein(*(at->resname[resnr])) ) {
731 /* mark this residue */
732 bResProcessed[resnr]=TRUE;
733 /* find out if this residue needs converting */
734 whatres=NOTSET;
735 for(j=0; j<resNR && whatres==NOTSET; j++) {
736 if ( ( !bPartial[j] &&
737 (strcasecmp(resnm,resnms[j])==0) ) ||
738 ( bPartial[j] &&
739 (strncasecmp(resnm,resnms[j],strlen(resnms[j]))==0) ) ) {
740 whatres=j;
741 /* get atoms we will be needing for the conversion */
742 nrfound=0;
743 for (k=0; atnms[j][k]; k++) {
744 ats[k]=NOTSET;
745 i0=i;
746 while (i<at->nr && at->atom[i].resnr==resnr && ats[k]==NOTSET) {
747 if (strcasecmp(*(at->atomname[i]),atnms[j][k])==0) {
748 ats[k]=i;
749 nrfound++;
750 }
751 i++;
752 }
753 /* if nothing found, search next atom from same point */
754 if (ats[k]==NOTSET)
755 i=i0;
756 }
757 /* now k is number of atom names in atnms[j] */
758 if (j==resHIS)
759 needed = k-3;
760 else
761 needed = k;
762 if (nrfound<needed)
763 fatal_error(0,"not enough atoms found (%d, need %d) in "
764 "residue %s %d while\n "
765 "generating aromatics dummy construction",
766 nrfound,needed,resnm,resnr+1);
767 }
768 }
769 /* the enums for every residue MUST correspond to atnms[residue] */
770 switch (whatres) {
771 case resPHE:
772 if (debug) fprintf(stderr,"PHE at %d\n",o2n[ats[0]]+1);
773 ndum+=gen_dums_phe(at, dummy_type, plist, nrfound, ats);
774 break;
775 case resTRP:
776 if (debug) fprintf(stderr,"TRP at %d\n",o2n[ats[0]]+1);
777 ndum+=gen_dums_trp(atype, &newx, &newatom, &newatomname, &o2n,
778 &newdummy_type, &newcgnr, symtab, &nadd, *x, cgnr,
779 at, dummy_type, plist, nrfound, ats, add_shift);
780 break;
781 case resTYR:
782 if (debug) fprintf(stderr,"TYR at %d\n",o2n[ats[0]]+1);
783 ndum+=gen_dums_tyr(at, dummy_type, plist, nrfound, ats);
784 break;
785 case resHIS:
786 if (debug) fprintf(stderr,"HIS at %d\n",o2n[ats[0]]+1);
787 ndum+=gen_dums_his(at, dummy_type, plist, nrfound, ats);
788 break;
789 case NOTSET:
790 /* this means this residue won't be processed */
791 break;
792 default:
793 fatal_error(0,"DEATH HORROR in do_dummies (%s:%d)",
794 __FILE__,__LINE__);
795 } /* switch whatres */
796 /* skip back to beginning of residue */
797 while(i>0 && at->atom[i-1].resnr==resnr)
798 i--;
799 } /* if bDummyAromatics & is protein */
800
801 /* now process the rest of the hydrogens */
802 /* only process hydrogen atoms which are not already set */
803 if ( ((*dummy_type)[i]==NOTSET) && is_hydrogen(*(at->atomname[i]))) {
804 /* find heavy atom, count #bonds from it and #H atoms bound to it
805 and return H atom numbers (Hatoms) and heavy atom numbers (heavies) */
806 count_bonds(i, &plist[F_BONDS], at->atomname,
807 &nrbonds, &nrHatoms, Hatoms, &Heavy, &nrheavies, heavies);
808 /* get Heavy atom type */
809 tpHeavy=get_atype(Heavy,at,nrtp,rtp);
810 strcpy(tpname,type2nm(tpHeavy,atype));
811 bWARNING=FALSE;
812 bAddDumParam=TRUE;
813 /* nested if's which check nrHatoms, nrbonds and tpname */
814 if (nrHatoms == 1) {
815 switch(nrbonds) {
816 case 2: /* -O-H */
817 (*dummy_type)[i]=F_BONDS;
818 break;
819 case 3: /* =CH-, -NH- or =NH+- */
820 (*dummy_type)[i]=F_DUMMY3FD;
821 break;
822 case 4: /* --CH- (tert) */
823 (*dummy_type)[i]=F_DUMMY4FD;
824 break;
825 default: /* nrbonds != 2, 3 or 4 */
826 bWARNING=TRUE;
827 }
828 } else if ( (nrHatoms == 2) && (nrbonds == 2) &&
829 (strcasecmp(tpname,"OW")==0) ) {
830 bAddDumParam=FALSE; /* this is water: skip these hydrogens */
831 if (bFirstWater) {
832 bFirstWater=FALSE;
833 if (debug)
834 fprintf(debug,
835 "Not converting hydrogens in water to dummy atoms\n");
836 }
837 } else if ( (nrHatoms == 2) && (nrbonds == 3) &&
838 (strcasecmp(tpname,"NL")!=0) ) {
839 /* =CH2 or =NH2 */
840 (*dummy_type)[Hatoms[0]] = F_DUMMY3FAD;
841 (*dummy_type)[Hatoms[1]] =-F_DUMMY3FAD;
842
843 } else if ( (nrHatoms == 2) && (nrbonds == 4) ) {
844 /* -CH2- */
845 (*dummy_type)[Hatoms[0]] = F_DUMMY3OUT;
846 (*dummy_type)[Hatoms[1]] =-F_DUMMY3OUT;
847
848 } else if ( ( (nrHatoms == 2) && (nrbonds == 3) &&
849 (strcasecmp(tpname,"NL")==0) ) ||
850 ( (nrHatoms == 3) && (nrbonds == 4) ) ) {
851 /* this tells how to set the hydrogen atoms */
852 int Hat_dummy_type[3] = { F_DUMMY3, F_DUMMY3OUT, F_DUMMY3OUT };
853 bool Hat_SwapParity[3] = { FALSE, TRUE, FALSE };
854
855 if (debug) fprintf(stderr,"-XH3 group at %d\n",i+1);
856 bAddDumParam=FALSE; /* we'll do this ourselves! */
857 /* -NH2 (umbrella), -NH3+ or -CH3 */
858 (*dummy_type)[Heavy] = F_DUMMY3;
859 for (j=0; j<nrHatoms; j++)
860 (*dummy_type)[Hatoms[j]] = Hat_dummy_type[j];
861 /* get dummy mass type from first char of heavy atom type (N or C) */
862 sprintf(name,"M%cH3",tpname[0]);
863 tpM=nm2type(name,atype);
864 /* make space for 2 masses: shift all atoms starting with 'Heavy' */
865 #define NMASS 2
866 i0=Heavy;
867 ni0=i0+nadd;
868 if (debug)
869 fprintf(stderr,"Inserting %d dummy masses at %d\n",NMASS,o2n[i0]+1);
870 nadd+=NMASS;
871 for(j=i0; j<at->nr; j++)
872 o2n[j]=j+nadd;
873 srenew(newx,at->nr+nadd);
874 srenew(newatom,at->nr+nadd);
875 srenew(newatomname,at->nr+nadd);
876 srenew(newdummy_type,at->nr+nadd);
877 srenew(newcgnr,at->nr+nadd);
878 for(j=0; j<NMASS; j++)
879 newatomname[at->nr+nadd-1-j]=NULL;
880
881 /* calculate starting position for the masses */
882 mHtot=0;
883 /* get atom masses, and set Heavy and Hatoms mass to zero */
884 for(j=0; j<nrHatoms; j++) {
885 mHtot += get_amass(Hatoms[j],at,nrtp,rtp);
886 at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
887 }
888 mtot = mHtot + get_amass(Heavy,at,nrtp,rtp);
889 at->atom[Heavy].m = at->atom[Heavy].mB = 0;
890 if (mHmult != 1.0)
891 mHtot *= mHmult;
892 fact2=mHtot/mtot;
893 fact=sqrt(fact2);
894 /* generate vectors parallel and perpendicular to rotational axis:
895 * rpar = Heavy -> Hcom
896 * rperp = Hcom -> H1 */
897 clear_rvec(rpar);
898 for(j=0; j<nrHatoms; j++)
899 rvec_inc(rpar,(*x)[Hatoms[j]]);
900 svmul(1.0/nrHatoms,rpar,rpar); /* rpar = ( H1+H2+H3 ) / 3 */
901 rvec_dec(rpar,(*x)[Heavy]); /* - Heavy */
902 rvec_sub((*x)[Hatoms[0]],(*x)[Heavy],rperp);
903 rvec_dec(rperp,rpar); /* rperp = H1 - Heavy - rpar */
904 /* calc mass positions */
905 svmul(fact2,rpar,temp);
906 for (j=0; (j<NMASS); j++) /* xM = xN + fact2 * rpar +/- fact * rperp */
907 rvec_add((*x)[Heavy],temp,newx[ni0+j]);
908 svmul(fact,rperp,temp);
909 rvec_inc(newx[ni0 ],temp);
910 rvec_dec(newx[ni0+1],temp);
911 /* set atom parameters for the masses */
912 for(j=0; (j<NMASS); j++) {
913 /* make name: "M??#" or "M?#" (? is atomname, # is number) */
914 name[0]='M';
915 for(k=0; (*at->atomname[Heavy])[k] && ( k < NMASS ); k++ )
916 name[k+1]=(*at->atomname[Heavy])[k];
917 name[k+1]=atomnamesuffix[j];
918 name[k+2]='\0';
919 newatomname[ni0+j] = put_symtab(symtab,name);
920 newatom[ni0+j].m = newatom[ni0+j].mB = mtot/NMASS;
921 newatom[ni0+j].q = newatom[ni0+j].qB = 0.0;
922 newatom[ni0+j].type = newatom[ni0+j].typeB = tpM;
923 newatom[ni0+j].ptype = eptAtom;
924 newatom[ni0+j].resnr = at->atom[i0].resnr;
925 newatom[ni0+j].chain = at->atom[i0].chain;
926 newdummy_type[ni0+j] = NOTSET;
927 newcgnr[ni0+j] = (*cgnr)[i0];
928 }
929 /* add constraints between dummy masses and to heavies[0] */
930 /* 'add_shift' says which atoms won't be renumbered afterwards */
931 my_add_param(&(plist[F_SHAKENC]),heavies[0], add_shift+ni0, NOTSET);
932 my_add_param(&(plist[F_SHAKENC]),heavies[0], add_shift+ni0+1,NOTSET);
933 my_add_param(&(plist[F_SHAKENC]),add_shift+ni0,add_shift+ni0+1,NOTSET);
934
935 /* generate Heavy, H1, H2 and H3 from M1, M2 and heavies[0] */
936 /* note that dummy_type cannot be NOTSET, because we just set it */
937 add_dum3_atoms (&plist[(*dummy_type)[Heavy]],
938 Heavy, heavies[0], add_shift+ni0, add_shift+ni0+1,
939 FALSE);
940 for(j=0; j<nrHatoms; j++)
941 add_dum3_atoms(&plist[(*dummy_type)[Hatoms[j]]],
942 Hatoms[j], heavies[0], add_shift+ni0, add_shift+ni0+1,
943 Hat_SwapParity[j]);
944 #undef NMASS
945 } else
946 bWARNING=TRUE;
947 if (bWARNING)
948 fprintf(stderr,
949 "Warning: cannot convert atom %d %s (bound to a heavy atom "
950 "type %s with \n"
951 " %d bonds and %d bound hydrogens atoms) to dummy "
952 "atom\n",
953 i+1,*(at->atomname[i]),tpname,nrbonds,nrHatoms);
954 if (bAddDumParam) {
955 /* add dummy parameters to topology,
956 also get rid of negative dummy_types */
957 add_dum_atoms(plist, (*dummy_type), Heavy, nrHatoms, Hatoms,
958 nrheavies, heavies);
959 /* transfer mass of dummy atom to Heavy atom */
960 for(j=0; j<nrHatoms; j++)
961 if (is_dum((*dummy_type)[Hatoms[j]])) {
962 at->atom[Heavy].m += at->atom[Hatoms[j]].m;
963 at->atom[Heavy].mB = at->atom[Heavy].m;
964 at->atom[Hatoms[j]].m = at->atom[Hatoms[j]].mB = 0;
965 }
966 }
967 ndum+=nrHatoms;
968 if (debug) {
969 fprintf(debug,"atom %d: ",o2n[i]+1);
970 print_bonds(debug,o2n,nrHatoms,Hatoms,Heavy,nrheavies,heavies);
971 }
972 } /* if dummy NOTSET & is hydrogen */
973
974 } /* for i < at->nr */
975
976 if (debug) {
977 fprintf(debug,"Before inserting new atoms:\n");
978 for(i=0; i<at->nr; i++)
979 fprintf(debug,"%4d %4d %4s %4d %4s %6d %-10s\n",i+1,o2n[i]+1,
980 at->atomname[i]?*(at->atomname[i]):"(NULL)",
981 at->atom[i].resnr,
982 at->resname[at->atom[i].resnr]?
983 *(at->resname[at->atom[i].resnr]):"(NULL)",
984 (*cgnr)[i],
985 ((*dummy_type)[i]==NOTSET) ?
986 "NOTSET" : interaction_function[(*dummy_type)[i]].name);
987 fprintf(debug,"new atoms to be inserted:\n");
988 for(i=0; i<at->nr+nadd; i++)
989 if (newatomname[i])
990 fprintf(debug,"%4d %4s %4d %6d %-10s\n",i+1,
991 newatomname[i]?*(newatomname[i]):"(NULL)",
992 newatom[i].resnr,newcgnr[i],
993 (newdummy_type[i]==NOTSET) ?
994 "NOTSET" : interaction_function[newdummy_type[i]].name);
995 }
996
997 /* add all original atoms to the new arrays, using o2n index array */
998 for(i=0; i<at->nr; i++) {
999 newatomname [o2n[i]] = at->atomname [i];
1000 newatom [o2n[i]] = at->atom [i];
1001 newdummy_type[o2n[i]] = (*dummy_type)[i];
1002 newcgnr [o2n[i]] = (*cgnr) [i];
1003 copy_rvec((*x)[i],newx[o2n[i]]);
1004 }
1005 /* throw away old atoms */
1006 sfree(at->atom);
1007 sfree(at->atomname);
1008 sfree(*dummy_type);
1009 sfree(*cgnr);
1010 sfree(*x);
1011 /* put in the new ones */
1012 at->nr += nadd;
1013 at->atom = newatom;
1014 at->atomname = newatomname;
1015 *dummy_type = newdummy_type;
1016 *cgnr = newcgnr;
1017 *x = newx;
1018 if (at->nr > add_shift)
1019 fatal_error(0,"Added impossible amount of dummy masses "
1020 "(%d on a total of %d atoms)\n",nadd,at->nr-nadd);
1021
1022 if (debug) {
1023 fprintf(debug,"After inserting new atoms:\n");
1024 for(i=0; i<at->nr; i++)
1025 fprintf(debug,"%4d %4s %4d %4s %6d %-10s\n",i+1,
1026 at->atomname[i]?*(at->atomname[i]):"(NULL)",
1027 at->atom[i].resnr,
1028 at->resname[at->atom[i].resnr]?
1029 *(at->resname[at->atom[i].resnr]):"(NULL)",
1030 (*cgnr)[i],
1031 ((*dummy_type)[i]==NOTSET) ?
1032 "NOTSET" : interaction_function[(*dummy_type)[i]].name);
1033 }
1034
1035 /* now renumber all the interactions because of the added atoms */
1036 for (ftype=0; ftype<F_NRE; ftype++) {
1037 params=&(plist[ftype]);
1038 if (debug)
1039 fprintf(debug,"Renumbering %d %s\n", params->nr,
1040 interaction_function[ftype].longname);
1041 for (i=0; i<params->nr; i++) {
1042 for (j=0; j<NRAL(ftype); j++)
1043 if (params->param[i].a[j]>=add_shift) {
1044 if (debug) fprintf(debug," [%u -> %u]",params->param[i].a[j],
1045 params->param[i].a[j]-add_shift);
1046 params->param[i].a[j]=params->param[i].a[j]-add_shift;
1047 } else {
1048 if (debug) fprintf(debug," [%u -> %d]",params->param[i].a[j],
1049 o2n[params->param[i].a[j]]);
1050 params->param[i].a[j]=o2n[params->param[i].a[j]];
1051 }
1052 if (debug) fprintf(debug,"\n");
1053 }
1054 }
1055 /* now check if atoms in the added constraints are in increasing order */
1056 params=&(plist[F_SHAKENC]);
1057 for(i=0; i<params->nr; i++)
1058 if ( params->param[i].AI > params->param[i].AJ ) {
1059 j = params->param[i].AJ;
1060 params->param[i].AJ = params->param[i].AI;
1061 params->param[i].AI = j;
1062 }
1063
1064 /* clean up */
1065 sfree(o2n);
1066
1067 /* tell the user what we did */
1068 fprintf(stderr,"Marked %d dummy atoms\n",ndum);
1069 fprintf(stderr,"Added %d dummy masses\n",nadd);
1070 fprintf(stderr,"Added %d new constraints\n",plist[F_SHAKENC].nr);
1071 }
1072
1073 void do_h_mass(t_params *psb, int dummy_type[], t_atoms *at, real mHmult)
1074 {
1075 int i,j,a;
1076
1077 /* loop over all atoms */
1078 for (i=0; i<at->nr; i++)
1079 /* adjust masses if i is hydrogen and not a dummy atom */
1080 if ( !is_dum(dummy_type[i]) && is_hydrogen(*(at->atomname[i])) ) {
1081 /* find bonded heavy atom */
1082 a=NOTSET;
1083 for(j=0; (j<psb->nr) && (a==NOTSET); j++) {
1084 /* if other atom is not a dummy, it is the one we want */
1085 if ( (psb->param[j].AI==i) &&
1086 !is_dum(dummy_type[psb->param[j].AJ]) )
1087 a=psb->param[j].AJ;
1088 else if ( (psb->param[j].AJ==i) &&
1089 !is_dum(dummy_type[psb->param[j].AI]) )
1090 a=psb->param[j].AI;
1091 }
1092 if (a==NOTSET)
1093 fatal_error(0,"Unbound hydrogen atom (%d) found while adjusting mass",
1094 i+1);
1095
1096 /* adjust mass of i (hydrogen) with mHmult
1097 and correct mass of a (bonded atom) with same amount */
1098 at->atom[a].m -= (mHmult-1.0)*at->atom[i].m;
1099 at->atom[a].mB-= (mHmult-1.0)*at->atom[i].m;
1100 at->atom[i].m *= mHmult;
1101 at->atom[i].mB*= mHmult;
1102 }
1103 }
1104
got merge a long time ago