?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
added Verlet scheme and NxN non-bonded functionality
[gromacs.git] / src / tools / dlist.c
blob7ff0660a54f44f0e1bd30a1375271ed9685f73d1
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 * Green Red Orange Magenta Azure Cyan Skyblue
34 */
35 #ifdef HAVE_CONFIG_H
36 #include <config.h>
37 #endif
38
39 #include <stdlib.h>
40
41 #include "string2.h"
42 #include "smalloc.h"
43 #include "gstat.h"
44 #include "gmx_fatal.h"
45 #include "index.h"
46
47 t_dlist *mk_dlist(FILE *log,
48 t_atoms *atoms, int *nlist,
49 gmx_bool bPhi, gmx_bool bPsi, gmx_bool bChi, gmx_bool bHChi,
50 int maxchi, int r0, gmx_residuetype_t rt)
51 {
52 int ires,i,j,k,ii;
53 t_dihatms atm,prev;
54 int nl=0,nc[edMax];
55 char *thisres;
56 t_dlist *dl;
57
58 snew(dl,atoms->nres+1);
59 prev.C = prev.O = -1;
60 for(i=0; (i<edMax); i++)
61 nc[i]=0;
62 ires = -1;
63 i = 0;
64 while (i<atoms->nr) {
65 ires = atoms->atom[i].resind;
66
67 /* Initiate all atom numbers to -1 */
68 atm.minC=atm.H=atm.N=atm.C=atm.O=atm.minO=-1;
69 for(j=0; (j<MAXCHI+3); j++)
70 atm.Cn[j]=-1;
71
72 /* Look for atoms in this residue */
73 /* maybe should allow for chis to hydrogens? */
74 while ((i<atoms->nr) && (atoms->atom[i].resind == ires)) {
75 if ((strcmp(*(atoms->atomname[i]),"H") == 0) ||
76 (strcmp(*(atoms->atomname[i]),"H1") == 0) )
77 atm.H=i;
78 else if (strcmp(*(atoms->atomname[i]),"N") == 0)
79 atm.N=i;
80 else if (strcmp(*(atoms->atomname[i]),"C") == 0)
81 atm.C=i;
82 else if ((strcmp(*(atoms->atomname[i]),"O") == 0) ||
83 (strcmp(*(atoms->atomname[i]),"O1") == 0))
84 atm.O=i;
85 else if (strcmp(*(atoms->atomname[i]),"CA") == 0)
86 atm.Cn[1]=i;
87 else if (strcmp(*(atoms->atomname[i]),"CB") == 0)
88 atm.Cn[2]=i;
89 else if ((strcmp(*(atoms->atomname[i]),"CG") == 0) ||
90 (strcmp(*(atoms->atomname[i]),"CG1") == 0) ||
91 (strcmp(*(atoms->atomname[i]),"OG") == 0) ||
92 (strcmp(*(atoms->atomname[i]),"OG1") == 0) ||
93 (strcmp(*(atoms->atomname[i]),"SG") == 0))
94 atm.Cn[3]=i;
95 else if ((strcmp(*(atoms->atomname[i]),"CD") == 0) ||
96 (strcmp(*(atoms->atomname[i]),"CD1") == 0) ||
97 (strcmp(*(atoms->atomname[i]),"SD") == 0) ||
98 (strcmp(*(atoms->atomname[i]),"OD1") == 0) ||
99 (strcmp(*(atoms->atomname[i]),"ND1") == 0))
100 atm.Cn[4]=i;
101 /* by grs - split the Cn[4] into 2 bits to check allowing dih to H */
102 else if (bHChi && ((strcmp(*(atoms->atomname[i]),"HG") == 0) ||
103 (strcmp(*(atoms->atomname[i]),"HG1") == 0)) )
104 atm.Cn[4]=i;
105 else if ((strcmp(*(atoms->atomname[i]),"CE") == 0) ||
106 (strcmp(*(atoms->atomname[i]),"CE1") == 0) ||
107 (strcmp(*(atoms->atomname[i]),"OE1") == 0) ||
108 (strcmp(*(atoms->atomname[i]),"NE") == 0))
109 atm.Cn[5]=i;
110 else if ((strcmp(*(atoms->atomname[i]),"CZ") == 0) ||
111 (strcmp(*(atoms->atomname[i]),"NZ") == 0))
112 atm.Cn[6]=i;
113 /* HChi flag here too */
114 else if (bHChi && (strcmp(*(atoms->atomname[i]),"NH1") == 0))
115 atm.Cn[7]=i;
116 i++;
117 }
118
119 thisres = *(atoms->resinfo[ires].name);
120
121 /* added by grs - special case for aromatics, whose chis above 2 are
122 not real and produce rubbish output - so set back to -1 */
123 if (strcmp(thisres,"PHE") == 0 ||
124 strcmp(thisres,"TYR") == 0 ||
125 strcmp(thisres,"PTR") == 0 ||
126 strcmp(thisres,"TRP") == 0 ||
127 strcmp(thisres,"HIS") == 0 ||
128 strcmp(thisres,"HISA") == 0 ||
129 strcmp(thisres,"HISB") == 0 ) {
130 for (ii=5 ; ii<=7 ; ii++)
131 atm.Cn[ii]=-1;
132 }
133 /* end fixing aromatics */
134
135 /* Special case for Pro, has no H */
136 if (strcmp(thisres,"PRO") == 0)
137 atm.H=atm.Cn[4];
138 /* Carbon from previous residue */
139 if (prev.C != -1)
140 atm.minC = prev.C;
141 if (prev.O != -1)
142 atm.minO = prev.O;
143 prev = atm;
144
145 /* Check how many dihedrals we have */
146 if ((atm.N != -1) && (atm.Cn[1] != -1) && (atm.C != -1) &&
147 (atm.O != -1) && ((atm.H != -1) || (atm.minC != -1))) {
148 dl[nl].resnr = ires+1;
149 dl[nl].atm = atm;
150 dl[nl].atm.Cn[0] = atm.N;
151 if ((atm.Cn[3] != -1) && (atm.Cn[2] != -1) && (atm.Cn[1] != -1)) {
152 nc[0]++;
153 if (atm.Cn[4] != -1) {
154 nc[1]++;
155 if (atm.Cn[5] != -1) {
156 nc[2]++;
157 if (atm.Cn[6] != -1) {
158 nc[3]++;
159 if (atm.Cn[7] != -1) {
160 nc[4]++;
161 if (atm.Cn[8] != -1) {
162 nc[5]++;
163 }
164 }
165 }
166 }
167 }
168 }
169 if ((atm.minC != -1) && (atm.minO != -1))
170 nc[6]++;
171 dl[nl].index=gmx_residuetype_get_index(rt,thisres);
172
173 sprintf(dl[nl].name,"%s%d",thisres,ires+r0);
174 nl++;
175 }
176 else if (debug)
177 fprintf(debug,"Could not find N atom but could find other atoms"
178 " in residue %s%d\n",thisres,ires+r0);
179 }
180 fprintf(stderr,"\n");
181 fprintf(log,"\n");
182 fprintf(log,"There are %d residues with dihedrals\n",nl);
183 j=0;
184 if (bPhi) j+=nl;
185 if (bPsi) j+=nl;
186 if (bChi)
187 for(i=0; (i<maxchi); i++)
188 j+=nc[i];
189 fprintf(log,"There are %d dihedrals\n",j);
190 fprintf(log,"Dihedral: ");
191 if (bPhi)
192 fprintf(log," Phi ");
193 if (bPsi)
194 fprintf(log," Psi ");
195 if (bChi)
196 for(i=0; (i<maxchi); i++)
197 fprintf(log,"Chi%d ",i+1);
198 fprintf(log,"\nNumber: ");
199 if (bPhi)
200 fprintf(log,"%4d ",nl);
201 if (bPsi)
202 fprintf(log,"%4d ",nl);
203 if (bChi)
204 for(i=0; (i<maxchi); i++)
205 fprintf(log,"%4d ",nc[i]);
206 fprintf(log,"\n");
207
208 *nlist=nl;
209
210 return dl;
211 }
212
213 gmx_bool has_dihedral(int Dih,t_dlist *dl)
214 {
215 gmx_bool b = FALSE;
216 int ddd;
217
218 switch (Dih) {
219 case edPhi:
220 b = ((dl->atm.H!=-1) && (dl->atm.N!=-1) && (dl->atm.Cn[1]!=-1) && (dl->atm.C!=-1));
221 break;
222 case edPsi:
223 b = ((dl->atm.N!=-1) && (dl->atm.Cn[1]!=-1) && (dl->atm.C!=-1) && (dl->atm.O!=-1));
224 break;
225 case edOmega:
226 b = ((dl->atm.minO!=-1) && (dl->atm.minC!=-1) && (dl->atm.N!=-1) && (dl->atm.Cn[1]!=-1));
227 break;
228 case edChi1:
229 case edChi2:
230 case edChi3:
231 case edChi4:
232 case edChi5:
233 case edChi6:
234 ddd = Dih - edChi1;
235 b = ((dl->atm.Cn[ddd]!=-1) && (dl->atm.Cn[ddd+1]!=-1)&&
236 (dl->atm.Cn[ddd+2]!=-1) && (dl->atm.Cn[ddd+3]!=-1));
237 break;
238 default:
239 pr_dlist(stdout,1,dl,1,0,TRUE,TRUE,TRUE,TRUE,MAXCHI);
240 gmx_fatal(FARGS,"Non existant dihedral %d in file %s, line %d",
241 Dih,__FILE__,__LINE__);
242 }
243 return b;
244 }
245
246 static void pr_one_ro(FILE *fp,t_dlist *dl,int nDih,real dt)
247 {
248 int k ;
249 for(k=0;k<NROT;k++)
250 fprintf(fp," %6.2f",dl->rot_occ[nDih][k]);
251 fprintf(fp,"\n");
252 }
253
254 static void pr_ntr_s2(FILE *fp,t_dlist *dl,int nDih,real dt)
255 {
256 fprintf(fp," %6.2f %6.2f\n",(dt == 0) ? 0 : dl->ntr[nDih]/dt,dl->S2[nDih]);
257 }
258
259 void pr_dlist(FILE *fp,int nl,t_dlist dl[],real dt, int printtype,
260 gmx_bool bPhi, gmx_bool bPsi,gmx_bool bChi,gmx_bool bOmega, int maxchi)
261 {
262 int i, Xi;
263
264 void (*pr_props)(FILE *, t_dlist *, int, real);
265
266 /* Analysis of dihedral transitions etc */
267
268 if (printtype == edPrintST){
269 pr_props=pr_ntr_s2 ;
270 fprintf(stderr,"Now printing out transitions and OPs...\n");
271 }else{
272 pr_props=pr_one_ro ;
273 fprintf(stderr,"Now printing out rotamer occupancies...\n");
274 fprintf(fp,"\nXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXXX\n\n");
275 }
276
277 /* change atom numbers from 0 based to 1 based */
278 for(i=0; (i<nl); i++) {
279 fprintf(fp,"Residue %s\n",dl[i].name);
280 if (printtype == edPrintST){
281 fprintf(fp," Angle [ AI, AJ, AK, AL] #tr/ns S^2D \n"
282 "--------------------------------------------\n");
283 } else {
284 fprintf(fp," Angle [ AI, AJ, AK, AL] rotamers 0 g(-) t g(+)\n"
285 "--------------------------------------------\n");
286 }
287 if (bPhi) {
288 fprintf(fp," Phi [%5d,%5d,%5d,%5d]",
289 (dl[i].atm.H == -1) ? 1+dl[i].atm.minC : 1+dl[i].atm.H,
290 1+dl[i].atm.N, 1+dl[i].atm.Cn[1], 1+dl[i].atm.C);
291 pr_props(fp,&dl[i],edPhi,dt);
292 }
293 if (bPsi) {
294 fprintf(fp," Psi [%5d,%5d,%5d,%5d]",1+dl[i].atm.N, 1+dl[i].atm.Cn[1],
295 1+dl[i].atm.C, 1+dl[i].atm.O);
296 pr_props(fp,&dl[i],edPsi,dt);
297 }
298 if (bOmega && has_dihedral(edOmega,&(dl[i]))) {
299 fprintf(fp," Omega [%5d,%5d,%5d,%5d]",1+dl[i].atm.minO, 1+dl[i].atm.minC,
300 1+dl[i].atm.N, 1+dl[i].atm.Cn[1]);
301 pr_props(fp,&dl[i],edOmega,dt);
302 }
303 for(Xi=0; Xi<MAXCHI; Xi++)
304 if (bChi && (Xi < maxchi) && (dl[i].atm.Cn[Xi+3] != -1) ) {
305 fprintf(fp," Chi%d[%5d,%5d,%5d,%5d]",Xi+1, 1+dl[i].atm.Cn[Xi],
306 1+dl[i].atm.Cn[Xi+1], 1+dl[i].atm.Cn[Xi+2],
307 1+dl[i].atm.Cn[Xi+3]);
308 pr_props(fp,&dl[i],Xi+edChi1,dt); /* Xi+2 was wrong here */
309 }
310 fprintf(fp,"\n");
311 }
312 }
313
314
315
316 int pr_trans(FILE *fp,int nl,t_dlist dl[],real dt,int Xi)
317 {
318 /* never called at the moment */
319
320 int i,nn,nz;
321
322 nz=0;
323 fprintf(fp,"\\begin{table}[h]\n");
324 fprintf(fp,"\\caption{Number of dihedral transitions per nanosecond}\n");
325 fprintf(fp,"\\begin{tabular}{|l|l|}\n");
326 fprintf(fp,"\\hline\n");
327 fprintf(fp,"Residue\t&$\\chi_%d$\t\\\\\n",Xi+1);
328 for(i=0; (i<nl); i++) {
329 nn=dl[i].ntr[Xi]/dt;
330
331 if (nn == 0) {
332 fprintf(fp,"%s\t&\\HL{%d}\t\\\\\n",dl[i].name,nn);
333 nz++;
334 }
335 else if (nn > 0)
336 fprintf(fp,"%s\t&\\%d\t\\\\\n",dl[i].name,nn);
337 }
338 fprintf(fp,"\\hline\n");
339 fprintf(fp,"\\end{tabular}\n");
340 fprintf(fp,"\\end{table}\n\n");
341
342 return nz;
343 }
344
The ultimate molecular dynamics simulation package