?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
changed reading hint
[gromacs/adressmacs.git] / src / tools / dtools.c
blobbf1863ced7105e550c1cb3e6459f085de9b52374
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_dtools_c = "$Id$";
30
31 #include "smalloc.h"
32 #include "strdb.h"
33 #include "futil.h"
34 #include "pdbio.h"
35 #include "vec.h"
36 #include "names.h"
37 #include "disco.h"
38 #include "pp2shift.h"
39
40 char *edc_names[edcNR+1] = { "NONBND", "BOND", "DISRE", NULL };
41
42 bool newres(int i,t_atom atom[])
43 {
44 return ((i == 0) || (atom[i].resnr != atom[i-1].resnr));
45 }
46
47 t_correct *init_corr(int maxnit,int nstprint,int nbcheck,int nstranlist,
48 int ngrow,bool bExplicit,bool bChiral,bool bPep,
49 bool bDump,real lowdev,bool bLowerOnly)
50 {
51 t_correct *c;
52
53 snew(c,1);
54 c->maxnit = maxnit;
55 c->nstprint = nstprint;
56 c->nbcheck = nbcheck;
57 c->nstranlist = nstranlist;
58 c->bExplicit = bExplicit;
59 c->bChiral = bChiral;
60 c->bPep = bPep;
61 c->bDump = bDump;
62 c->lodev = lowdev;
63 c->maxdist = 0;
64 c->ndist = 0;
65 c->ngrow = ngrow;
66 c->bLowerOnly = bLowerOnly;
67
68 return c;
69 }
70
71 void make_tags(t_correct *c,int natom)
72 {
73 int i,n,ai;
74
75 if (c->ndist == 0)
76 fatal_error(0,"Can not make tags without distances. %s, %d",
77 __FILE__,__LINE__);
78
79 /* First sort the distances */
80
81 /* Now make the tags */
82 snew(c->tag,natom+1);
83 ai = c->d[0].ai;
84 n = 0;
85 for(i=0; (i<c->ndist); i++) {
86 if (c->d[i].ai != ai) {
87 /* New tag starts here, but skip over possible missing atoms */
88 while ((n < natom) && (n<c->d[i].ai)) {
89 n++;
90 c->tag[n] = i;
91 ai = c->d[i].ai;
92 }
93 /* else
94 fatal_error(0,"Too many atoms, or distances not sorted");*/
95 }
96 }
97 if (n < natom) {
98 while (n < natom) {
99 n++;
100 c->tag[n] = i;
101 }
102 }
103 else
104 fatal_error(0,"Too many atoms, or distances not sorted");
105 fprintf(stderr,"There are %d tags for %d atoms\n",n,natom);
106
107 if (debug)
108 for(i=0; (i<=natom); i++)
109 fprintf(debug,"tag[%5d] = %d\n",i,c->tag[i]);
110 }
111
112 void center_in_box(int natom,rvec xin[],matrix box,rvec xout[])
113 {
114 int i,m;
115 rvec xcm,dx;
116
117 /* Dump the optimization trajectory to an xtc file */
118 /* Put the whole thing in the center of the box 1st */
119 clear_rvec(xcm);
120 for(i=0; (i<natom); i++) {
121 copy_rvec(xin[i],xout[i]);
122 rvec_inc(xcm,xin[i]);
123 }
124 for(m=0; (m<DIM); m++)
125 dx[m] = 0.5*box[m][m]-xcm[m]/natom;
126 for(i=0; (i<natom); i++)
127 rvec_inc(xout[i],dx);
128 }
129
130 void define_peptide_bonds(FILE *log,t_atoms *atoms,t_correct *c)
131 {
132 int i,npep,naa,nlist;
133 char **aa;
134 t_dlist *dlist;
135
136 naa = get_strings("aminoacids.dat",&aa);
137 dlist = mk_dlist(log,atoms,&nlist,TRUE,TRUE,FALSE,0,1,naa,aa);
138 for(i=0; (i<naa); i++)
139 sfree(aa[i]);
140 sfree(aa);
141
142 npep = 0;
143 snew(c->pepbond,nlist);
144 snew(c->omega,nlist);
145 for(i=0; (i<nlist); i++) {
146 if (has_dihedral(edOmega,&dlist[i])) {
147 c->pepbond[npep].ai = dlist[i].atm.minC;
148 c->pepbond[npep].aj = dlist[i].atm.minO;
149 c->pepbond[npep].ak = dlist[i].atm.N;
150 c->pepbond[npep].al = dlist[i].atm.H;
151 c->pepbond[npep].am = dlist[i].atm.Cn[1];
152 npep++;
153 }
154 }
155 c->npep = npep;
156 if (debug)
157 pr_dlist(debug,nlist,dlist,1.0);
158 sfree(dlist);
159 fprintf(log,"There are %d peptide bonds\n",npep);
160 }
161
162 static char *aname(t_atoms *atoms,int i)
163 {
164 static char buf[32];
165
166 sprintf(buf,"%4s%3d:%4s",
167 *(atoms->resname[atoms->atom[i].resnr]),
168 atoms->atom[i].resnr+1,
169 *(atoms->atomname[i]));
170
171 return buf;
172 }
173
174 int find_atom(t_atoms *atoms,int resnr,int j0,char *nm)
175 {
176 int j,aa=-1;
177
178 for(j=j0; (j < atoms->nr) && (atoms->atom[j].resnr == resnr); j++)
179 if (strcmp(*atoms->atomname[j],nm) == 0) {
180 aa = j;
181 break;
182 }
183 return aa;
184 }
185
186 void define_impropers(FILE *log,t_atoms *atoms,t_correct *c)
187 {
188 typedef struct {
189 char *res,*aa[4];
190 } t_impdef;
191
192 t_impdef id[] = {
193 { NULL, "CA", "N", "C", "CB" },
194 { NULL, "N", "CA", "H1", "H3" },
195 { "LYSH", "NZ", "CE", "HZ2", "HZ1" },
196 { "LYSH", "NZ", "CE", "HZ1", "HZ3" },
197 { "LEU", "CG", "CD2", "CD1", "CB" },
198 { "VAL", "CB", "CG2", "CG1", "CA" },
199 { "ILE", "CB", "CG1", "CG2", "CA" },
200 { "THR", "CB", "OG1", "CG2", "CA" }
201 };
202 #define NID asize(id)
203 int i,j,k,l,aa[4],nimp;
204
205 nimp = 0;
206 for(i=0; (i<atoms->nres); i++) {
207 /* Skip until residue number */
208 for(j=0; (j<atoms->nr) && (atoms->atom[j].resnr != i); j++)
209 ;
210 for(k=0; (k<NID); k++) {
211 if ((id[k].res == NULL) ||
212 (strcmp(id[k].res,*atoms->resname[i]) == 0)) {
213 /* This (i) is the right residue to look for this improper (k) */
214 for(l=0; (l<4); l++)
215 aa[l] = find_atom(atoms,i,j,id[k].aa[l]);
216 if ((aa[0] != -1) && (aa[1] != -1) && (aa[2] != -1) && (aa[3] != -1)) {
217 srenew(c->imp,nimp+1);
218 srenew(c->idih,nimp+1);
219 c->imp[nimp].ai = aa[0];
220 c->imp[nimp].aj = aa[1];
221 c->imp[nimp].ak = aa[2];
222 c->imp[nimp].al = aa[3];
223 nimp++;
224 }
225 }
226 }
227 }
228 c->nimp = nimp;
229
230 fprintf(log,"There are %d impropers\n",c->nimp);
231 if (debug) {
232 fprintf(debug,"Overview of improper dihedrals\n");
233 for(i=0; (i<c->nimp); i++) {
234 fprintf(debug," %s",aname(atoms,c->imp[i].ai));
235 fprintf(debug," %s",aname(atoms,c->imp[i].aj));
236 fprintf(debug," %s",aname(atoms,c->imp[i].ak));
237 fprintf(debug," %s\n",aname(atoms,c->imp[i].al));
238 }
239 }
240 }
241
242 void pr_dist(FILE *fp,bool bHeader,t_correct *c,int i)
243 {
244 real ideal=0;
245
246 if (bHeader)
247 fprintf(fp,"#%4s%5s%10s%10s%10s\n","ai","aj","ideal","lb","ub");
248 switch (c->d[i].cons_type) {
249 case edcBOND:
250 ideal = 5*(c->d[i].lb+c->d[i].ub);
251 break;
252 case edcNONBOND:
253 ideal = 0;
254 break;
255 case edcDISRE:
256 ideal = -1;
257 break;
258 default:
259 fatal_error(0,"cons_type for distance %d = %d\n",i,c->d[i].cons_type);
260 }
261 fprintf(fp,"%5d%5d%10.5f%10.5f%10.5f\n",1+c->d[i].ai,1+c->d[i].aj,
262 ideal,10*c->d[i].lb,10*c->d[i].ub);
263 }
264
265 void pr_distances(FILE *fp,t_correct *c)
266 {
267 int i;
268
269 for(i=0; (i<c->ndist); i++)
270 pr_dist(fp,(i==0),c,i);
271 }
272
273 void add_dist(t_correct *c,int ai,int aj,real lb,real ideal,real ub,real w[])
274 {
275 int n = c->ndist;
276
277 if ((w[ai] != 0) || (w[aj] != 0)) {
278 if (n == c->maxdist) {
279 c->maxdist += 100;
280 srenew(c->d,c->maxdist);
281 srenew(c->ip,c->maxdist);
282 }
283 c->d[n].ai = ai;
284 c->d[n].aj = aj;
285 if (ideal > 0)
286 c->d[n].cons_type = edcBOND;
287 else if (ideal == 0.0)
288 c->d[n].cons_type = edcNONBOND;
289 else
290 c->d[n].cons_type = edcDISRE;
291 c->d[n].lb = lb;
292 c->d[n].ub = ub;
293 c->d[n].wi = w[ai]/(w[ai]+w[aj]);
294 c->ip[n] = n;
295 c->ndist++;
296 }
297 }
298
299 void define_dist(FILE *log,t_topology *top,t_correct *c,real weight[])
300 {
301 int i,k,nra,ndist;
302 int ai,aj,tp,ftype;
303 real lb=0,ub=0,b0;
304
305 c->ndist = 0;
306 c->d = NULL;
307 c->ip = NULL;
308 snew(c->bViol,top->atoms.nr);
309 for(ftype=0; (ftype<F_NRE); ftype++) {
310 if ((interaction_function[ftype].flags & IF_BOND) ||
311 (ftype == F_DISRES)) {
312 nra = interaction_function[ftype].nratoms+1;
313
314 ndist = top->idef.il[ftype].nr/nra;
315
316 for(i=k=0; (i<ndist); i++,k+=nra) {
317 tp = top->idef.il[ftype].iatoms[k];
318 ai = top->idef.il[ftype].iatoms[k+1];
319 aj = top->idef.il[ftype].iatoms[k+2];
320 b0 = 0;
321 switch (ftype) {
322 case F_DISRES:
323 lb = top->idef.iparams[tp].disres.low;
324 ub = top->idef.iparams[tp].disres.up1;
325 break;
326 case F_SHAKE:
327 b0 = top->idef.iparams[tp].shake.dA;
328 break;
329 case F_BONDS:
330 case F_G96BONDS:
331 b0 = top->idef.iparams[tp].harmonic.rA;
332 break;
333 case F_MORSE:
334 b0 = top->idef.iparams[tp].morse.b0;
335 break;
336 }
337 if (b0 != 0.0) {
338 lb = b0*0.98;
339 ub = b0*1.02;
340 }
341 /* CHECK THIS! */
342 b0 = (lb + ub)*0.5;
343 add_dist(c,ai,aj,lb,b0,ub,weight);
344 }
345 }
346 }
347 }
348
349 void read_dist(FILE *log,char *fn,int natom,t_correct *c,real weight[])
350 {
351 FILE *fp;
352 char buf[256];
353 int ai,aj,i,nline=0;
354 double ideal,lb,ub;
355 int nnn[edcNR];
356
357 for(i=0; (i<edcNR); i++)
358 nnn[i] = 0;
359
360 snew(c->bViol,natom);
361
362 fp = ffopen(fn,"r");
363 while (fgets(buf,255,fp)) {
364 nline ++;
365 if (buf[0] != '#') {
366 if (sscanf(buf,"%d%d%lf%lf%lf",&ai,&aj,&ideal,&lb,&ub) != 5)
367 fprintf(stderr,"Error in %s, line %d\n",fn,nline);
368 else {
369 add_dist(c,ai-1,aj-1,lb*0.1,ideal*0.1,ub*0.1,weight);
370 nnn[c->d[c->ndist-1].cons_type]++;
371 }
372 }
373 }
374 ffclose(fp);
375
376 fprintf(log,"Read distances from file %s\n",fn);
377 for(i=0; (i<edcNR); i++)
378 fprintf(log,"There were %d %s distances\n",
379 nnn[i],edc_names[i]);
380 }
381
382 real *read_weights(char *fn,int natom)
383 {
384 t_atoms newatoms;
385 int i;
386 rvec *x;
387 matrix box;
388 char title[256];
389 real *w;
390
391 /* Read the weights from the occupancy field in the pdb file */
392 snew(x,natom);
393 init_t_atoms(&newatoms,natom,TRUE);
394 read_pdb_conf(fn,title,&newatoms,x,box,FALSE);
395 snew(w,newatoms.nr);
396 for(i=0; (i<newatoms.nr); i++)
397 w[i] = newatoms.pdbinfo[i].occup;
398 free_t_atoms(&newatoms);
399 sfree(x);
400
401 return w;
402 }
403
404 void pr_corr(FILE *log,t_correct *c)
405 {
406 fprintf(log,"Parameters for this disco run\n");
407 fprintf(log,"maxnit = %d\n",c->maxnit);
408 fprintf(log,"nbcheck = %d\n",c->nbcheck);
409 fprintf(log,"nstprint = %d\n",c->nstprint);
410 fprintf(log,"nstranlist = %d\n",c->nstranlist);
411 fprintf(log,"ngrow = %d\n",c->ngrow);
412 fprintf(log,"bExplicit = %s\n",BOOL(c->bExplicit));
413 fprintf(log,"bChiral = %s\n",BOOL(c->bChiral));
414 fprintf(log,"bPep = %s\n",BOOL(c->bPep));
415 fprintf(log,"bDump = %s\n",BOOL(c->bDump));
416 fprintf(log,"ndist = %d\n",c->ndist);
417 fprintf(log,"npep = %d\n",c->npep);
418 fprintf(log,"nimp = %d\n",c->nimp);
419 fprintf(log,"lowdev = %g\n",c->lodev);
420 fflush(log);
421 }
422
423 void measure_dist(FILE *log,int natom,rvec x[],t_correct *c,real weight[],
424 real cutoff)
425 {
426 int ai,aj;
427 rvec dx;
428 real ideal;
429
430 snew(c->bViol,natom);
431
432 for(ai=0; (ai<natom); ai++)
433 for(aj=ai+1; (aj<natom); aj++) {
434 rvec_sub(x[ai],x[aj],dx);
435 ideal = 10.0*norm(dx);
436 if ((ideal < cutoff) || (cutoff == 0))
437 add_dist(c,ai,aj,ideal*0.98,ideal,ideal*1.02,weight);
438 }
439 }
440
441 void check_dist(FILE *log,t_correct *c)
442 {
443 int i;
444 real tol=0.001;
445
446 fprintf(log,"Checking distances for internal consistency\n");
447 for(i=0; (i<c->ndist); i++) {
448 if ((c->d[i].ub - c->d[i].lb) < tol) {
449 pr_dist(log,TRUE,c,i);
450 }
451 }
452 }
453
got merge a long time ago