?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / kernel / pdb2top.c
blob88938e0ff731b6a91bda6b0cc1db1dd39a572bab
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_pdb2top_c = "$Id$";
30
31 #include <stdio.h>
32 #include <math.h>
33 #include "vec.h"
34 #include "copyrite.h"
35 #include "assert.h"
36 #include "smalloc.h"
37 #include "macros.h"
38 #include "symtab.h"
39 #include "futil.h"
40 #include "fatal.h"
41 #include "pdb2top.h"
42 #include "topexcl.h"
43 #include "topdirs.h"
44 #include "toputil.h"
45 #include "h_db.h"
46 #include "pgutil.h"
47 #include "resall.h"
48 #include "topio.h"
49 #include "physics.h"
50 #include "pdbio.h"
51 #include "gen_ad.h"
52 #include "filenm.h"
53 #include "index.h"
54 #include "gen_dum.h"
55 #include "add_par.h"
56
57 /* this must correspond to enum in pdb2top.h */
58 char *hh[ehisNR] = { "HISA", "HISB", "HISH", "HIS1" };
59
60 static bool missing_atoms(t_restp *rp, int resnr,
61 t_atoms *at, int i0, int i, bool bCTer)
62 {
63 int j,k;
64 char *name;
65 bool bFound, bRet;
66
67 bRet=FALSE;
68 for (j=0; j<rp->natom; j++) {
69 name=*(rp->atomname[j]);
70 if ((name[0]!='H') && (name[0]!='h') && (!bCTer || (name[0]!='O'))) {
71 bFound=FALSE;
72 for (k=i0; k<i; k++)
73 bFound=(bFound || !strcasecmp(*(at->atomname[k]),name));
74 if (!bFound) {
75 bRet=TRUE;
76 fprintf(stderr,"\nWARNING: "
77 "atom %s is missing in residue %s %d in the pdb file\n\n",
78 name,*(at->resname[resnr]),resnr+1);
79 }
80 }
81 }
82
83 return bRet;
84 }
85
86 bool is_int(double x)
87 {
88 const double tol = 1e-4;
89 int ix;
90
91 if (x < 0)
92 x=-x;
93 ix=gmx_nint(x);
94
95 return (fabs(x-ix) < tol);
96 }
97
98 char *choose_ff(void)
99 {
100 typedef struct { char *desc,*fn; } t_fff;
101 static char *fnsel;
102 FILE *in;
103 t_fff *fff;
104 int i,nff,sel;
105 char *c,buf[STRLEN],fn[32];
106
107 in=libopen("FF.dat");
108 fgets2(buf,255,in);
109 sscanf(buf,"%d",&nff);
110 snew(fff,nff);
111 for(i=0; (i<nff); i++) {
112 fgets2(buf,255,in);
113 sscanf(buf,"%s",fn);
114 fff[i].fn=strdup(fn);
115 /* Search for next non-space character, there starts description */
116 c=&(buf[strlen(fn)+1]);
117 while (isspace(*c)) c++;
118 fff[i].desc=strdup(c);
119 }
120 fclose(in);
121
122 if (nff > 1) {
123 printf("\nSelect the Force Field:\n");
124 for(i=0; (i<nff); i++)
125 printf("%2d: %s\n",i,fff[i].desc);
126 do {
127 fgets(buf,STRLEN,stdin);
128 sscanf(buf,"%d",&sel);
129 } while ((sel < 0) || (sel >= nff));
130 }
131 else
132 sel=0;
133
134 fnsel=strdup(fff[sel].fn);
135
136 for(i=0; (i<nff); i++) {
137 sfree(fff[i].desc);
138 sfree(fff[i].fn);
139 }
140 sfree(fff);
141
142 return fnsel;
143 }
144
145 static void name2type(t_atoms *at, int **cgnr, t_atomtype *atype,
146 t_restp restp[])
147 {
148 int i,j,prevresnr,resnr,i0,prevcg,cg,curcg;
149 char *name;
150 bool bProt, bNterm;
151 double qt;
152
153 resnr=-1;
154 bProt=FALSE;
155 bNterm=FALSE;
156 i0=0;
157 snew(*cgnr,at->nr);
158 qt=0;
159 prevcg=NOTSET;
160 curcg=0;
161 cg=-1;
162 j=NOTSET;
163 for(i=0; (i<at->nr); i++) {
164 prevresnr=resnr;
165 if (at->atom[i].resnr != resnr) {
166 resnr=at->atom[i].resnr;
167 bProt=is_protein(*(at->resname[resnr]));
168 bNterm=bProt && (resnr == 0);
169 if (resnr>0)
170 missing_atoms(&restp[prevresnr],resnr,at,i0,i,
171 (!bProt && is_protein(restp[prevresnr].resname)));
172 i0=i;
173 }
174 if (at->atom[i].m == 0) {
175 if (debug)
176 fprintf(debug,"atom %d%s: curcg=%d, prevcg=%d, cg=%d\n",
177 i+1,*(at->atomname[i]),curcg,prevcg,
178 j==NOTSET ? NOTSET : restp[resnr].cgnr[j]);
179 qt=0;
180 prevcg=cg;
181 name=*(at->atomname[i]);
182 j=search_jtype(&restp[resnr],name,bNterm);
183 at->atom[i].type = restp[resnr].atom[j].type;
184 at->atom[i].q = restp[resnr].atom[j].q;
185 at->atom[i].m = atype->atom[restp[resnr].atom[j].type].m;
186 cg = restp[resnr].cgnr[j];
187 if ( (cg != prevcg) || (resnr != prevresnr) )
188 curcg++;
189 } else {
190 if (debug)
191 fprintf(debug,"atom %d%s: curcg=%d, qt=%g, is_int=%d\n",
192 i+1,*(at->atomname[i]),curcg,qt,is_int(qt));
193 cg=-1;
194 if (is_int(qt)) {
195 qt=0;
196 curcg++;
197 }
198 qt+=at->atom[i].q;
199 }
200 (*cgnr)[i]=curcg;
201 at->atom[i].typeB = at->atom[i].type;
202 at->atom[i].qB = at->atom[i].q;
203 at->atom[i].mB = at->atom[i].m;
204 }
205 missing_atoms(&restp[resnr],resnr,at,i0,i,
206 (!bProt || is_protein(restp[resnr].resname)));
207 }
208
209 static void print_top_heavy_H(FILE *out, real mHmult)
210 {
211 if (mHmult!=1.0) {
212 fprintf(out,"; heavy hydrogens:\n");
213 if (mHmult==4.0)
214 fprintf(out,"#define HEAVY_H\n\n");
215 else
216 fprintf(stderr,"WARNING: unsupported proton mass multiplier (%g) "
217 "in pdb2top\n",mHmult);
218 }
219 }
220
221 void print_top_comment(FILE *out, char *title, bool bITP)
222 {
223 fprintf(out,"; This is your %stopology file\n",bITP ? "include " : "");
224 fprintf(out,"; %s\n\n",title[0]?title:cool_quote());
225 }
226
227 void print_top_header(FILE *out, char *title, bool bITP, char *ff, real mHmult)
228 {
229 print_top_comment(out,title,bITP);
230
231 print_top_heavy_H(out, mHmult);
232 fprintf(out,"; Include forcefield parameters\n");
233 fprintf(out,"#include \"%s.itp\"\n\n",ff);
234 }
235
236 static void print_top_posre(FILE *out,char *pr)
237 {
238 fprintf(out,"; Include Position restraint file\n");
239 fprintf(out,"#ifdef POSRES\n");
240 fprintf(out,"#include \"%s\"\n",pr);
241 fprintf(out,"#endif\n\n");
242 }
243
244 static void print_top_water(FILE *out)
245 {
246 fprintf(out,"; Include water topology\n");
247 fprintf(out,"#ifdef FLEX_SPC\n");
248 fprintf(out,"#include \"flexspc.itp\"\n");
249 fprintf(out,"#else\n");
250 fprintf(out,"#include \"spc.itp\"\n");
251 fprintf(out,"#endif\n");
252 fprintf(out,"\n");
253 fprintf(out,"#ifdef POSRES_WATER\n");
254 fprintf(out,"; Position restraint for each water oxygen\n");
255 fprintf(out,"[ position_restraints ]\n");
256 fprintf(out,";%3s %5s %9s %10s %10s\n","i","funct","fcx","fcy","fcz");
257 fprintf(out,"%4d %4d %10g %10g %10g\n",1,1,1000.0,1000.0,1000.0);
258 fprintf(out,"#endif\n");
259 fprintf(out,"\n");
260 }
261
262 static void print_top_system(FILE *out, char *title)
263 {
264 fprintf(out,"[ %s ]\n",dir2str(d_system));
265 fprintf(out,"; Name\n");
266 fprintf(out,"%s\n\n",title[0]?title:"Protein");
267 }
268
269 void print_top_mols(FILE *out, char *title,
270 int nincl, char **incls, int nmol, t_mols *mols)
271 {
272 int i;
273
274 if (nincl>0) {
275 fprintf(out,"; Include chain topologies\n");
276 for (i=0; (i<nincl); i++)
277 fprintf(out,"#include \"%s\"\n",incls[i]);
278 fprintf(out,"\n");
279 }
280
281 print_top_water(out);
282 print_top_system(out, title);
283
284 if (nmol) {
285 fprintf(out,"[ %s ]\n",dir2str(d_molecules));
286 fprintf(out,"; %-15s %5s\n","Compound","#mols");
287 for (i=0; (i<nmol); i++)
288 fprintf(out,"%-15s %5d\n",mols[i].name,mols[i].nr);
289 }
290 }
291
292 void write_top(FILE *out, char *pr,char *molname,
293 t_atoms *at,int bts[],t_params plist[],t_block *excl,
294 t_atomtype *atype,int *cgnr, int nrexcl)
295 /* NOTE: nrexcl is not the size of *excl! */
296 {
297 if (at && atype && cgnr) {
298 fprintf(out,"[ %s ]\n",dir2str(d_moleculetype));
299 fprintf(out,"; %-15s %5s\n","Name","nrexcl");
300 fprintf(out,"%-15s %5d\n\n",molname?molname:"Protein",nrexcl);
301
302 print_atoms(out, atype, at, cgnr);
303 print_bondeds(out,at->nr,d_bonds, F_BONDS, bts[ebtsBONDS], plist);
304 print_bondeds(out,at->nr,d_constraints,F_SHAKE, 0, plist);
305 print_bondeds(out,at->nr,d_constraints,F_SHAKENC, 0, plist);
306 print_bondeds(out,at->nr,d_pairs, F_LJ14, 0, plist);
307 print_bondeds(out,at->nr,d_angles, F_ANGLES, bts[ebtsANGLES],plist);
308 print_bondeds(out,at->nr,d_dihedrals, F_PDIHS, bts[ebtsPDIHS], plist);
309 print_bondeds(out,at->nr,d_dihedrals, F_IDIHS, bts[ebtsIDIHS], plist);
310 print_bondeds(out,at->nr,d_dum3, F_DUMMY3, 0, plist);
311 print_bondeds(out,at->nr,d_dum3, F_DUMMY3FD, 0, plist);
312 print_bondeds(out,at->nr,d_dum3, F_DUMMY3FAD,0, plist);
313 print_bondeds(out,at->nr,d_dum3, F_DUMMY3OUT,0, plist);
314 print_bondeds(out,at->nr,d_dum4, F_DUMMY4FD, 0, plist);
315
316 if ( excl && excl->nr > 0 )
317 print_excl(out,excl);
318
319 if (pr)
320 print_top_posre(out,pr);
321 }
322 }
323
324 static atom_id search_res_atom(char *type,int resnr,
325 int natom,t_atom at[],char **aname[])
326 {
327 int i;
328
329 for(i=0; (i<natom); i++)
330 if (at[i].resnr == resnr)
331 return search_atom(type,i,natom,at,aname);
332
333 return NO_ATID;
334 }
335
336 static void do_ssbonds(t_params *ps,int natoms,t_atom atom[],char **aname[],
337 int nssbonds,t_ssbond *ssbonds)
338 {
339 int i,ri,rj;
340 atom_id ai,aj;
341
342 for(i=0; (i<nssbonds); i++) {
343 ri = ssbonds[i].res1;
344 rj = ssbonds[i].res2;
345 ai = search_res_atom(ssbonds[i].a1,ri,natoms,atom,aname);
346 aj = search_res_atom(ssbonds[i].a2,rj,natoms,atom,aname);
347 if ((ai == NO_ATID) || (aj == NO_ATID))
348 fatal_error(0,"Trying to make impossible special bond (%s-%s)!",
349 ssbonds[i].a1,ssbonds[i].a2);
350 add_param(ps,ai,aj,NULL,NULL);
351 }
352 }
353
354 static void at2bonds(t_params *psb, t_hackblock *hb,
355 int natoms, t_atom atom[], char **aname[],
356 int nres, rvec x[],
357 real long_bond_dist, real short_bond_dist)
358 {
359 int resnr,i,j,k;
360 atom_id ai,aj;
361 real dist2, long_bond_dist2, short_bond_dist2;
362
363 long_bond_dist2 = sqr(long_bond_dist);
364 short_bond_dist2 = sqr(short_bond_dist);
365
366 fprintf(stderr,"Making bonds...\n");
367 i=0;
368 for(resnr=0; (resnr < nres) && (i<natoms); resnr++) {
369 /* add bonds from list of bonded interactions */
370 for(j=0; j < hb[resnr].rb[ebtsBONDS].nb; j++) {
371 ai=search_atom(hb[resnr].rb[ebtsBONDS].b[j].AI,i,natoms,atom,aname);
372 aj=search_atom(hb[resnr].rb[ebtsBONDS].b[j].AJ,i,natoms,atom,aname);
373 if ( ai != NO_ATID && aj != NO_ATID ) {
374 dist2 = distance2(x[ai],x[aj]);
375 if (dist2 > long_bond_dist2 )
376 fprintf(stderr,"Warning: Long Bond (%d-%d = %g nm)\n",
377 ai+1,aj+1,sqrt(dist2));
378 else if (dist2 < short_bond_dist2 )
379 fprintf(stderr,"Warning: Short Bond (%d-%d = %g nm)\n",
380 ai+1,aj+1,sqrt(dist2));
381 add_param(psb,ai,aj,NULL,hb[resnr].rb[ebtsBONDS].b[j].s);
382 }
383 }
384 /* add bonds from list of hacks (each added atom gets a bond) */
385 while( (i<natoms) && (atom[i].resnr == resnr) ) {
386 for(j=0; j < hb[resnr].nhack; j++)
387 if ( ( hb[resnr].hack[j].tp > 0 ||
388 hb[resnr].hack[j].oname==NULL ) &&
389 strcmp(hb[resnr].hack[j].AI,*(aname[i])) == 0 ) {
390 switch(hb[resnr].hack[j].tp) {
391 case 9: /* COOH terminus */
392 add_param(psb,i,i+1,NULL,NULL); /* C-O */
393 add_param(psb,i,i+2,NULL,NULL); /* C-OA */
394 add_param(psb,i+2,i+3,NULL,NULL); /* OA-H */
395 break;
396 default:
397 for(k=0; (k<hb[resnr].hack[j].nr); k++)
398 add_param(psb,i,i+k+1,NULL,NULL);
399 }
400 }
401 i++;
402 }
403 /* we're now at the start of the next residue */
404 }
405 }
406
407 static int pcompar(const void *a, const void *b)
408 {
409 t_param *pa,*pb;
410 int d;
411 pa=(t_param *)a;
412 pb=(t_param *)b;
413
414 d = pa->AI - pb->AI;
415 if (d == 0) {
416 d = pa->AJ - pb->AJ;
417 if (d == 0)
418 /* we'll keep the first bond in the list,
419 doing inverse sort will put the bond with the longest string first */
420 d = -strcmp(pa->s, pb->s);
421 }
422 return d;
423 }
424
425 static void clean_bonds(t_params *ps)
426 {
427 int i,j;
428 atom_id a;
429
430 if (ps->nr > 0) {
431 /* swap atomnumbers in bond if first larger than second: */
432 for(i=0; (i<ps->nr); i++)
433 if ( ps->param[i].AJ < ps->param[i].AI ) {
434 a = ps->param[i].AI;
435 ps->param[i].AI = ps->param[i].AJ;
436 ps->param[i].AJ = a;
437 }
438
439 /* Sort bonds */
440 qsort(ps->param,ps->nr,(size_t)sizeof(ps->param[0]),pcompar);
441
442 /* remove doubles */
443 j=1;
444 for (i=1; (i<ps->nr); i++) {
445 if ( ps->param[i].AI != ps->param[j-1].AI ||
446 ps->param[i].AJ != ps->param[j-1].AJ ) {
447 ps->param[j] = ps->param[i];
448 j++;
449 } else
450 sfree(ps->param[i].s);
451 }
452 fprintf(stderr,"Number of bonds was %d, now %d\n",ps->nr,j);
453 ps->nr=j;
454 } else
455 fprintf(stderr,"No bonds\n");
456 }
457
458 void print_sums(t_atoms *atoms, bool bSystem)
459 {
460 double m,qtot;
461 int i;
462 char *where;
463
464 if (bSystem)
465 where=" in system";
466 else
467 where="";
468
469 m=0;
470 qtot=0;
471 for(i=0; (i<atoms->nr); i++) {
472 m+=atoms->atom[i].m;
473 qtot+=atoms->atom[i].q;
474 }
475
476 fprintf(stderr,"Total mass%s %.3f a.m.u.\n",where,m);
477 fprintf(stderr,"Total charge%s %.3f e\n",where,qtot);
478 }
479
480 void get_hackblocks_rtp(t_hackblock **hb, t_restp **restp,
481 int nrtp, t_restp rtp[], int nres, char **resname[],
482 t_hackblock *ntdb, t_hackblock *ctdb, int rn, int rc)
483 {
484 int i, j, k, l;
485 t_restp *res;
486 char buf[STRLEN];
487 char Hnum[6]="123456";
488
489 snew(*hb,nres);
490 snew(*restp,nres);
491 /* first the termini */
492 if (rn>=0)
493 copy_t_hackblock(ntdb, &(*hb)[rn]);
494 if (rc>=0)
495 merge_t_hackblock(ctdb, &(*hb)[rc]);
496
497 /* then the whole rtp */
498 for(i=0; i < nres; i++) {
499 res = search_rtp(*resname[i],nrtp,rtp);
500 copy_t_restp(res, &(*restp)[i]);
501 merge_t_bondeds(res->rb, (*hb)[i].rb);
502 }
503
504 /* now perform t_hack's on t_restp's,
505 i.e. add's and deletes from termini database will be
506 added to/removed from residue topology
507 we'll do this on one big dirty loop, so it won't make easy reading! */
508 for(i=0; i < nres; i++)
509 for(j=0; j < (*hb)[i].nhack; j++)
510 if ( (*hb)[i].hack[j].nr ) {
511 /* find atom in restp */
512 for(l=0; l < (*restp)[i].natom; l++)
513 if ( ( (*hb)[i].hack[j].oname==NULL &&
514 strcmp((*hb)[i].hack[j].AI, *(*restp)[i].atomname[l])==0 ) ||
515 ( (*hb)[i].hack[j].oname!=NULL &&
516 strcmp((*hb)[i].hack[j].oname,*(*restp)[i].atomname[l])==0 ) )
517 break;
518 if (l == (*restp)[i].natom)
519 fatal_error(0,"atom %s not found in residue %d%s "
520 "while combining tdb and rtp",
521 (*hb)[i].hack[j].oname!=NULL ?
522 (*hb)[i].hack[j].oname : (*hb)[i].hack[j].AI,
523 i+1,*resname[i]);
524 if ( (*hb)[i].hack[j].oname==NULL ) {
525 /* we're adding: */
526 if (debug)
527 fprintf(debug,"adding atom(s) %s to atom %s in res %d%s in rtp\n",
528 (*hb)[i].hack[j].nname,
529 *(*restp)[i].atomname[l], i+1, (*restp)[i].resname);
530 strcpy(buf, (*hb)[i].hack[j].nname);
531 buf[strlen(buf)+1]='\0';
532 if ( (*hb)[i].hack[j].nr > 1 )
533 buf[strlen(buf)]='-';
534 /* make space */
535 (*restp)[i].natom += (*hb)[i].hack[j].nr;
536 srenew((*restp)[i].atom, (*restp)[i].natom);
537 srenew((*restp)[i].atomname, (*restp)[i].natom);
538 srenew((*restp)[i].cgnr, (*restp)[i].natom);
539 /* shift the rest */
540 for(k=(*restp)[i].natom-1; k > l+(*hb)[i].hack[j].nr; k--) {
541 (*restp)[i].atom[k] =
542 (*restp)[i].atom [k - (*hb)[i].hack[j].nr];
543 (*restp)[i].atomname[k] =
544 (*restp)[i].atomname[k - (*hb)[i].hack[j].nr];
545 (*restp)[i].cgnr[k] =
546 (*restp)[i].cgnr [k - (*hb)[i].hack[j].nr];
547 }
548 /* now add them */
549 for(k=0; k < (*hb)[i].hack[j].nr; k++) {
550 /* set counter in atomname */
551 if ( (*hb)[i].hack[j].nr > 1 )
552 buf[strlen(buf)-1]=Hnum[k];
553 snew( (*restp)[i].atomname[l+1+k], 1);
554 (*restp)[i].atom [l+1+k] = *(*hb)[i].hack[j].atom;
555 *(*restp)[i].atomname[l+1+k] = strdup(buf);
556 if ( (*hb)[i].hack[j].cgnr != NOTSET )
557 (*restp)[i].cgnr [l+1+k] = (*hb)[i].hack[j].cgnr;
558 else
559 (*restp)[i].cgnr [l+1+k] = (*restp)[i].cgnr[l];
560 }
561 } else /* oname != NULL */
562 if ( (*hb)[i].hack[j].nname==NULL ) {
563 /* we're deleting */
564 if (debug)
565 fprintf(debug, "deleting atom %s from res %d%s in rtp\n",
566 *(*restp)[i].atomname[l],
567 i+1,(*restp)[i].resname);
568 /* shift the rest */
569 (*restp)[i].natom--;
570 for(k=l; k < (*restp)[i].natom; k++) {
571 (*restp)[i].atom [k] = (*restp)[i].atom [k+1];
572 (*restp)[i].atomname[k] = (*restp)[i].atomname[k+1];
573 (*restp)[i].cgnr [k] = (*restp)[i].cgnr [k+1];
574 }
575 /* give back space */
576 srenew((*restp)[i].atom, (*restp)[i].natom);
577 srenew((*restp)[i].atomname, (*restp)[i].natom);
578 srenew((*restp)[i].cgnr, (*restp)[i].natom);
579 } else { /* nname != NULL */
580 /* we're replacing */
581 if (debug)
582 fprintf(debug, "replacing atom %s by %s in res %d%s in rtp\n",
583 *(*restp)[i].atomname[l], (*hb)[i].hack[j].nname,
584 i+1,(*restp)[i].resname);
585 snew( (*restp)[i].atomname[l], 1);
586 (*restp)[i].atom[l] = *(*hb)[i].hack[j].atom;
587 *(*restp)[i].atomname[l] = strdup((*hb)[i].hack[j].nname);
588 if ( (*hb)[i].hack[j].cgnr != NOTSET )
589 (*restp)[i].cgnr [l] = (*hb)[i].hack[j].cgnr;
590 }
591 }
592 }
593
594 void pdb2top(FILE *top_file, char *posre_fn, char *molname,
595 t_atoms *atoms, rvec **x, t_atomtype *atype, t_symtab *tab,
596 int bts[], int nrtp, t_restp rtp[],
597 t_hackblock *ntdb, t_hackblock *ctdb,
598 bool bH14, int rn, int rc, bool bAlldih,
599 bool bDummies, bool bDummyAromatics, real mHmult,
600 int nssbonds, t_ssbond *ssbonds, int nrexcl,
601 real long_bond_dist, real short_bond_dist)
602 {
603 t_hackblock *hb;
604 t_restp *restp;
605 t_params plist[F_NRE];
606 t_nextnb nnb;
607 int *cgnr;
608 int *dummy_type;
609 int i;
610
611 init_plist(plist);
612
613 /* lookup hackblocks and rtp for all residues */
614 get_hackblocks_rtp(&hb, &restp, nrtp, rtp, atoms->nres, atoms->resname,
615 ntdb, ctdb, rn, rc);
616 /* ideally, now we would not need the rtp itself anymore, but do
617 everything using the hb and restp arrays. Unfortunately, that
618 requires some re-thinking of code in gen_dum.c, which I won't
619 do now :( AF 26-7-99 */
620
621 if (debug) {
622 print_resall(debug, bts, atoms->nres, restp, atype);
623 dump_hb(debug, atoms->nres, hb);
624 }
625
626 /* Make bonds */
627 at2bonds(&(plist[F_BONDS]), hb,
628 atoms->nr, atoms->atom, atoms->atomname, atoms->nres, *x,
629 long_bond_dist, short_bond_dist);
630
631 /* specbonds: disulphide bonds & heme-his */
632 do_ssbonds(&(plist[F_BONDS]),
633 atoms->nr, atoms->atom, atoms->atomname, nssbonds, ssbonds);
634
635 name2type(atoms, &cgnr, atype, restp);
636
637 /* Cleanup bonds (sort and rm doubles) */
638 clean_bonds(&(plist[F_BONDS]));
639
640 snew(dummy_type,atoms->nr);
641 for(i=0; i<atoms->nr; i++)
642 dummy_type[i]=NOTSET;
643 if (bDummies)
644 /* determine which atoms will be dummies and add dummy masses
645 also renumber atom numbers in plist[0..F_NRE]! */
646 do_dummies(nrtp, rtp, atype, atoms, tab, x, plist,
647 &dummy_type, &cgnr, mHmult, bDummyAromatics);
648
649 /* Make Angles and Dihedrals */
650 fprintf(stderr,"Generating angles and dihedrals...\n");
651 init_nnb(&nnb,atoms->nr,4);
652 gen_nnb(&nnb,plist);
653 print_nnb(&nnb,"NNB");
654 gen_pad(&nnb,atoms,bH14,plist,hb,bAlldih);
655 done_nnb(&nnb);
656
657 /* set mass of all remaining hydrogen atoms */
658 if (mHmult != 1.0)
659 do_h_mass(&(plist[F_BONDS]),dummy_type,atoms,mHmult);
660 sfree(dummy_type);
661
662 /* Cleanup bonds (sort and rm doubles) */
663 clean_bonds(&(plist[F_BONDS]));
664
665 fprintf(stderr,
666 "There are %4d dihedrals, %4d impropers, %4d angles\n"
667 " %4d pairs, %4d bonds and %4d dummies\n",
668 plist[F_PDIHS].nr, plist[F_IDIHS].nr, plist[F_ANGLES].nr,
669 plist[F_LJ14].nr, plist[F_BONDS].nr,
670 plist[F_DUMMY2].nr +
671 plist[F_DUMMY3].nr +
672 plist[F_DUMMY3FD].nr +
673 plist[F_DUMMY3FAD].nr +
674 plist[F_DUMMY3OUT].nr +
675 plist[F_DUMMY4FD].nr );
676
677 print_sums(atoms, FALSE);
678
679 if (top_file) {
680 fprintf(stderr,"Writing topology\n");
681 write_top(top_file, posre_fn, molname,
682 atoms, bts, plist, NULL, atype, cgnr, nrexcl);
683 }
684
685 /* cleaning up */
686 free_t_hackblock(atoms->nres, &hb);
687 free_t_restp(atoms->nres, &restp);
688
689 /* we should clean up hb and restp here, but that is a *L*O*T* of work! */
690 sfree(cgnr);
691 for (i=0; (i<F_NRE); i++)
692 sfree(plist[i].param);
693 }
got merge a long time ago