2 * This file is part of the GROMACS molecular simulation package.
4 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
5 * Copyright (c) 2001-2004, The GROMACS development team.
6 * Copyright (c) 2013,2014, by the GROMACS development team, led by
7 * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl,
8 * and including many others, as listed in the AUTHORS file in the
9 * top-level source directory and at http://www.gromacs.org.
11 * GROMACS is free software; you can redistribute it and/or
12 * modify it under the terms of the GNU Lesser General Public License
13 * as published by the Free Software Foundation; either version 2.1
14 * of the License, or (at your option) any later version.
16 * GROMACS is distributed in the hope that it will be useful,
17 * but WITHOUT ANY WARRANTY; without even the implied warranty of
18 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
19 * Lesser General Public License for more details.
21 * You should have received a copy of the GNU Lesser General Public
22 * License along with GROMACS; if not, see
23 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
24 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
26 * If you want to redistribute modifications to GROMACS, please
27 * consider that scientific software is very special. Version
28 * control is crucial - bugs must be traceable. We will be happy to
29 * consider code for inclusion in the official distribution, but
30 * derived work must not be called official GROMACS. Details are found
31 * in the README & COPYING files - if they are missing, get the
32 * official version at http://www.gromacs.org.
34 * To help us fund GROMACS development, we humbly ask that you cite
35 * the research papers on the package. Check out http://www.gromacs.org.
41 #include "gromacs/utility/smalloc.h"
44 #include "gromacs/fileio/xvgr.h"
46 #include "gromacs/commandline/pargs.h"
47 #include "gromacs/math/vec.h"
48 #include "gromacs/topology/index.h"
49 #include "gromacs/fileio/pdbio.h"
50 #include "gromacs/fileio/tpxio.h"
51 #include "gromacs/fileio/trxio.h"
52 #include "gromacs/utility/futil.h"
54 #include "gromacs/pbcutil/rmpbc.h"
55 #include "gromacs/fileio/confio.h"
58 #include "gromacs/linearalgebra/eigensolver.h"
59 #include "gromacs/math/do_fit.h"
61 static real
find_pdb_bfac(t_atoms
*atoms
, t_resinfo
*ri
, char *atomnm
)
66 strcpy(rresnm
, *ri
->name
);
68 for (i
= 0; (i
< atoms
->nr
); i
++)
70 if ((ri
->nr
== atoms
->resinfo
[atoms
->atom
[i
].resind
].nr
) &&
71 (ri
->ic
== atoms
->resinfo
[atoms
->atom
[i
].resind
].ic
) &&
72 (strcmp(*atoms
->resinfo
[atoms
->atom
[i
].resind
].name
, rresnm
) == 0) &&
73 (strstr(*atoms
->atomname
[i
], atomnm
) != NULL
))
80 fprintf(stderr
, "\rCan not find %s%d-%s in pdbfile\n",
81 rresnm
, ri
->nr
, atomnm
);
85 return atoms
->pdbinfo
[i
].bfac
;
88 void correlate_aniso(const char *fn
, t_atoms
*ref
, t_atoms
*calc
,
89 const output_env_t oenv
)
94 fp
= xvgropen(fn
, "Correlation between X-Ray and Computed Uij", "X-Ray",
96 for (i
= 0; (i
< ref
->nr
); i
++)
98 if (ref
->pdbinfo
[i
].bAnisotropic
)
100 for (j
= U11
; (j
<= U23
); j
++)
102 fprintf(fp
, "%10d %10d\n", ref
->pdbinfo
[i
].uij
[j
], calc
->pdbinfo
[i
].uij
[j
]);
109 static void average_residues(double f
[], double **U
, int uind
,
110 int isize
, atom_id index
[], real w_rls
[],
119 for (i
= 0; i
< isize
; i
++)
121 av
+= w_rls
[index
[i
]]*(f
!= NULL
? f
[i
] : U
[i
][uind
]);
122 m
+= w_rls
[index
[i
]];
124 atoms
->atom
[index
[i
]].resind
!= atoms
->atom
[index
[i
+1]].resind
)
129 for (j
= start
; j
<= i
; j
++)
136 for (j
= start
; j
<= i
; j
++)
148 void print_dir(FILE *fp
, real
*Uaver
)
150 real eigvec
[DIM
*DIM
];
155 fprintf(fp
, "MSF X Y Z\n");
156 for (d
= 0; d
< DIM
; d
++)
158 fprintf(fp
, " %c ", 'X'+d
-XX
);
159 for (m
= 0; m
< DIM
; m
++)
161 fprintf(fp
, " %9.2e", Uaver
[3*m
+d
]);
163 fprintf(fp
, "%s\n", m
== DIM
? " (nm^2)" : "");
166 for (m
= 0; m
< DIM
*DIM
; m
++)
172 eigensolver(tmp
, DIM
, 0, DIM
, eigval
, eigvec
);
174 fprintf(fp
, "\n Eigenvectors\n\n");
175 fprintf(fp
, "Eigv %-8.2e %-8.2e %-8.2e (nm^2)\n\n",
176 eigval
[2], eigval
[1], eigval
[0]);
177 for (d
= 0; d
< DIM
; d
++)
179 fprintf(fp
, " %c ", 'X'+d
-XX
);
180 for (m
= DIM
-1; m
>= 0; m
--)
182 fprintf(fp
, "%7.4f ", eigvec
[3*m
+d
]);
188 int gmx_rmsf(int argc
, char *argv
[])
190 const char *desc
[] = {
191 "[THISMODULE] computes the root mean square fluctuation (RMSF, i.e. standard ",
192 "deviation) of atomic positions in the trajectory (supplied with [TT]-f[tt])",
193 "after (optionally) fitting to a reference frame (supplied with [TT]-s[tt]).[PAR]",
194 "With option [TT]-oq[tt] the RMSF values are converted to B-factor",
195 "values, which are written to a [TT].pdb[tt] file with the coordinates, of the",
196 "structure file, or of a [TT].pdb[tt] file when [TT]-q[tt] is specified.",
197 "Option [TT]-ox[tt] writes the B-factors to a file with the average",
199 "With the option [TT]-od[tt] the root mean square deviation with",
200 "respect to the reference structure is calculated.[PAR]",
201 "With the option [TT]-aniso[tt], [THISMODULE] will compute anisotropic",
202 "temperature factors and then it will also output average coordinates",
203 "and a [TT].pdb[tt] file with ANISOU records (corresonding to the [TT]-oq[tt]",
204 "or [TT]-ox[tt] option). Please note that the U values",
205 "are orientation-dependent, so before comparison with experimental data",
206 "you should verify that you fit to the experimental coordinates.[PAR]",
207 "When a [TT].pdb[tt] input file is passed to the program and the [TT]-aniso[tt]",
209 "a correlation plot of the Uij will be created, if any anisotropic",
210 "temperature factors are present in the [TT].pdb[tt] file.[PAR]",
211 "With option [TT]-dir[tt] the average MSF (3x3) matrix is diagonalized.",
212 "This shows the directions in which the atoms fluctuate the most and",
215 static gmx_bool bRes
= FALSE
, bAniso
= FALSE
, bdevX
= FALSE
, bFit
= TRUE
;
217 { "-res", FALSE
, etBOOL
, {&bRes
},
218 "Calculate averages for each residue" },
219 { "-aniso", FALSE
, etBOOL
, {&bAniso
},
220 "Compute anisotropic termperature factors" },
221 { "-fit", FALSE
, etBOOL
, {&bFit
},
222 "Do a least squares superposition before computing RMSF. Without this you must make sure that the reference structure and the trajectory match." }
225 int step
, nre
, natoms
, i
, g
, m
, teller
= 0;
226 real t
, lambda
, *w_rls
, *w_rms
;
232 t_atoms
*pdbatoms
, *refatoms
;
236 rvec
*x
, *pdbx
, *xref
;
243 FILE *fp
; /* the graphics file */
244 const char *devfn
, *dirfn
;
252 real bfac
, pdb_bfac
, *Uaver
;
256 double *rmsf
, invcount
, totmass
;
260 gmx_rmpbc_t gpbc
= NULL
;
264 const char *leg
[2] = { "MD", "X-Ray" };
267 { efTRX
, "-f", NULL
, ffREAD
},
268 { efTPS
, NULL
, NULL
, ffREAD
},
269 { efNDX
, NULL
, NULL
, ffOPTRD
},
270 { efPDB
, "-q", NULL
, ffOPTRD
},
271 { efPDB
, "-oq", "bfac", ffOPTWR
},
272 { efPDB
, "-ox", "xaver", ffOPTWR
},
273 { efXVG
, "-o", "rmsf", ffWRITE
},
274 { efXVG
, "-od", "rmsdev", ffOPTWR
},
275 { efXVG
, "-oc", "correl", ffOPTWR
},
276 { efLOG
, "-dir", "rmsf", ffOPTWR
}
278 #define NFILE asize(fnm)
280 if (!parse_common_args(&argc
, argv
, PCA_CAN_TIME
| PCA_CAN_VIEW
| PCA_BE_NICE
,
281 NFILE
, fnm
, asize(pargs
), pargs
, asize(desc
), desc
, 0, NULL
,
287 bReadPDB
= ftp2bSet(efPDB
, NFILE
, fnm
);
288 devfn
= opt2fn_null("-od", NFILE
, fnm
);
289 dirfn
= opt2fn_null("-dir", NFILE
, fnm
);
291 read_tps_conf(ftp2fn(efTPS
, NFILE
, fnm
), title
, &top
, &ePBC
, &xref
, NULL
, box
, TRUE
);
292 snew(w_rls
, top
.atoms
.nr
);
294 fprintf(stderr
, "Select group(s) for root mean square calculation\n");
295 get_index(&top
.atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
), 1, &isize
, &index
, &grpnames
);
298 for (i
= 0; i
< isize
; i
++)
300 w_rls
[index
[i
]] = top
.atoms
.atom
[index
[i
]].m
;
303 /* Malloc the rmsf arrays */
304 snew(xav
, isize
*DIM
);
306 for (i
= 0; i
< isize
; i
++)
318 get_stx_coordnum(opt2fn("-q", NFILE
, fnm
), &npdbatoms
);
321 init_t_atoms(pdbatoms
, npdbatoms
, TRUE
);
322 init_t_atoms(refatoms
, npdbatoms
, TRUE
);
323 snew(pdbx
, npdbatoms
);
324 /* Read coordinates twice */
325 read_stx_conf(opt2fn("-q", NFILE
, fnm
), title
, pdbatoms
, pdbx
, NULL
, NULL
, pdbbox
);
326 read_stx_conf(opt2fn("-q", NFILE
, fnm
), title
, refatoms
, pdbx
, NULL
, NULL
, pdbbox
);
330 pdbatoms
= &top
.atoms
;
331 refatoms
= &top
.atoms
;
333 npdbatoms
= pdbatoms
->nr
;
334 snew(pdbatoms
->pdbinfo
, npdbatoms
);
335 copy_mat(box
, pdbbox
);
340 sub_xcm(xref
, isize
, index
, top
.atoms
.atom
, xcm
, FALSE
);
343 natom
= read_first_x(oenv
, &status
, ftp2fn(efTRX
, NFILE
, fnm
), &t
, &x
, box
);
347 gpbc
= gmx_rmpbc_init(&top
.idef
, ePBC
, natom
);
350 /* Now read the trj again to compute fluctuations */
356 /* Remove periodic boundary */
357 gmx_rmpbc(gpbc
, natom
, box
, x
);
359 /* Set center of mass to zero */
360 sub_xcm(x
, isize
, index
, top
.atoms
.atom
, xcm
, FALSE
);
362 /* Fit to reference structure */
363 do_fit(natom
, w_rls
, xref
, x
);
366 /* Calculate Anisotropic U Tensor */
367 for (i
= 0; i
< isize
; i
++)
370 for (d
= 0; d
< DIM
; d
++)
372 xav
[i
*DIM
+ d
] += x
[aid
][d
];
373 for (m
= 0; m
< DIM
; m
++)
375 U
[i
][d
*DIM
+ m
] += x
[aid
][d
]*x
[aid
][m
];
382 /* Calculate RMS Deviation */
383 for (i
= 0; (i
< isize
); i
++)
386 for (d
= 0; (d
< DIM
); d
++)
388 rmsd_x
[i
][d
] += sqr(x
[aid
][d
]-xref
[aid
][d
]);
395 while (read_next_x(oenv
, status
, &t
, x
, box
));
400 gmx_rmpbc_done(gpbc
);
404 invcount
= 1.0/count
;
405 snew(Uaver
, DIM
*DIM
);
407 for (i
= 0; i
< isize
; i
++)
409 for (d
= 0; d
< DIM
; d
++)
411 xav
[i
*DIM
+ d
] *= invcount
;
413 for (d
= 0; d
< DIM
; d
++)
415 for (m
= 0; m
< DIM
; m
++)
417 U
[i
][d
*DIM
+ m
] = U
[i
][d
*DIM
+ m
]*invcount
418 - xav
[i
*DIM
+ d
]*xav
[i
*DIM
+ m
];
419 Uaver
[3*d
+m
] += top
.atoms
.atom
[index
[i
]].m
*U
[i
][d
*DIM
+ m
];
422 totmass
+= top
.atoms
.atom
[index
[i
]].m
;
424 for (d
= 0; d
< DIM
*DIM
; d
++)
431 for (d
= 0; d
< DIM
*DIM
; d
++)
433 average_residues(NULL
, U
, d
, isize
, index
, w_rls
, &top
.atoms
);
439 for (i
= 0; i
< isize
; i
++)
442 pdbatoms
->pdbinfo
[aid
].bAnisotropic
= TRUE
;
443 pdbatoms
->pdbinfo
[aid
].uij
[U11
] = 1e6
*U
[i
][XX
*DIM
+ XX
];
444 pdbatoms
->pdbinfo
[aid
].uij
[U22
] = 1e6
*U
[i
][YY
*DIM
+ YY
];
445 pdbatoms
->pdbinfo
[aid
].uij
[U33
] = 1e6
*U
[i
][ZZ
*DIM
+ ZZ
];
446 pdbatoms
->pdbinfo
[aid
].uij
[U12
] = 1e6
*U
[i
][XX
*DIM
+ YY
];
447 pdbatoms
->pdbinfo
[aid
].uij
[U13
] = 1e6
*U
[i
][XX
*DIM
+ ZZ
];
448 pdbatoms
->pdbinfo
[aid
].uij
[U23
] = 1e6
*U
[i
][YY
*DIM
+ ZZ
];
460 for (i
= 0; i
< isize
; i
++)
462 rmsf
[i
] = U
[i
][XX
*DIM
+ XX
] + U
[i
][YY
*DIM
+ YY
] + U
[i
][ZZ
*DIM
+ ZZ
];
467 fprintf(stdout
, "\n");
468 print_dir(stdout
, Uaver
);
469 fp
= gmx_ffopen(dirfn
, "w");
470 print_dir(fp
, Uaver
);
474 for (i
= 0; i
< isize
; i
++)
480 /* Write RMSF output */
483 bfac
= 8.0*M_PI
*M_PI
/3.0*100;
484 fp
= xvgropen(ftp2fn(efXVG
, NFILE
, fnm
), "B-Factors",
485 label
, "(A\\b\\S\\So\\N\\S2\\N)", oenv
);
486 xvgr_legend(fp
, 2, leg
, oenv
);
487 for (i
= 0; (i
< isize
); i
++)
489 if (!bRes
|| i
+1 == isize
||
490 top
.atoms
.atom
[index
[i
]].resind
!= top
.atoms
.atom
[index
[i
+1]].resind
)
492 resind
= top
.atoms
.atom
[index
[i
]].resind
;
493 pdb_bfac
= find_pdb_bfac(pdbatoms
, &top
.atoms
.resinfo
[resind
],
494 *(top
.atoms
.atomname
[index
[i
]]));
496 fprintf(fp
, "%5d %10.5f %10.5f\n",
497 bRes
? top
.atoms
.resinfo
[top
.atoms
.atom
[index
[i
]].resind
].nr
: index
[i
]+1, rmsf
[i
]*bfac
,
505 fp
= xvgropen(ftp2fn(efXVG
, NFILE
, fnm
), "RMS fluctuation", label
, "(nm)", oenv
);
506 for (i
= 0; i
< isize
; i
++)
508 if (!bRes
|| i
+1 == isize
||
509 top
.atoms
.atom
[index
[i
]].resind
!= top
.atoms
.atom
[index
[i
+1]].resind
)
511 fprintf(fp
, "%5d %8.4f\n",
512 bRes
? top
.atoms
.resinfo
[top
.atoms
.atom
[index
[i
]].resind
].nr
: index
[i
]+1, sqrt(rmsf
[i
]));
518 for (i
= 0; i
< isize
; i
++)
520 pdbatoms
->pdbinfo
[index
[i
]].bfac
= 800*M_PI
*M_PI
/3.0*rmsf
[i
];
525 for (i
= 0; i
< isize
; i
++)
527 rmsf
[i
] = (rmsd_x
[i
][XX
]+rmsd_x
[i
][YY
]+rmsd_x
[i
][ZZ
])/count
;
531 average_residues(rmsf
, NULL
, 0, isize
, index
, w_rls
, &top
.atoms
);
533 /* Write RMSD output */
534 fp
= xvgropen(devfn
, "RMS Deviation", label
, "(nm)", oenv
);
535 for (i
= 0; i
< isize
; i
++)
537 if (!bRes
|| i
+1 == isize
||
538 top
.atoms
.atom
[index
[i
]].resind
!= top
.atoms
.atom
[index
[i
+1]].resind
)
540 fprintf(fp
, "%5d %8.4f\n",
541 bRes
? top
.atoms
.resinfo
[top
.atoms
.atom
[index
[i
]].resind
].nr
: index
[i
]+1, sqrt(rmsf
[i
]));
547 if (opt2bSet("-oq", NFILE
, fnm
))
549 /* Write a .pdb file with B-factors and optionally anisou records */
550 for (i
= 0; i
< isize
; i
++)
552 rvec_inc(xref
[index
[i
]], xcm
);
554 write_sto_conf_indexed(opt2fn("-oq", NFILE
, fnm
), title
, pdbatoms
, pdbx
,
555 NULL
, ePBC
, pdbbox
, isize
, index
);
557 if (opt2bSet("-ox", NFILE
, fnm
))
559 /* Misuse xref as a temporary array */
560 for (i
= 0; i
< isize
; i
++)
562 for (d
= 0; d
< DIM
; d
++)
564 xref
[index
[i
]][d
] = xcm
[d
] + xav
[i
*DIM
+ d
];
567 /* Write a .pdb file with B-factors and optionally anisou records */
568 write_sto_conf_indexed(opt2fn("-ox", NFILE
, fnm
), title
, pdbatoms
, xref
, NULL
,
569 ePBC
, pdbbox
, isize
, index
);
573 correlate_aniso(opt2fn("-oc", NFILE
, fnm
), refatoms
, pdbatoms
, oenv
);
574 do_view(oenv
, opt2fn("-oc", NFILE
, fnm
), "-nxy");
576 do_view(oenv
, opt2fn("-o", NFILE
, fnm
), "-nxy");
579 do_view(oenv
, opt2fn("-od", NFILE
, fnm
), "-nxy");