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,2015,2016,2017, 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.
43 #include "gromacs/commandline/pargs.h"
44 #include "gromacs/commandline/viewit.h"
45 #include "gromacs/fileio/confio.h"
46 #include "gromacs/fileio/pdbio.h"
47 #include "gromacs/fileio/trxio.h"
48 #include "gromacs/fileio/xvgr.h"
49 #include "gromacs/gmxana/gmx_ana.h"
50 #include "gromacs/gmxana/princ.h"
51 #include "gromacs/linearalgebra/eigensolver.h"
52 #include "gromacs/math/do_fit.h"
53 #include "gromacs/math/functions.h"
54 #include "gromacs/math/utilities.h"
55 #include "gromacs/math/vec.h"
56 #include "gromacs/pbcutil/rmpbc.h"
57 #include "gromacs/topology/index.h"
58 #include "gromacs/topology/topology.h"
59 #include "gromacs/utility/arraysize.h"
60 #include "gromacs/utility/cstringutil.h"
61 #include "gromacs/utility/futil.h"
62 #include "gromacs/utility/smalloc.h"
64 static real
find_pdb_bfac(const t_atoms
*atoms
, t_resinfo
*ri
, char *atomnm
)
69 std::strcpy(rresnm
, *ri
->name
);
71 for (i
= 0; (i
< atoms
->nr
); i
++)
73 if ((ri
->nr
== atoms
->resinfo
[atoms
->atom
[i
].resind
].nr
) &&
74 (ri
->ic
== atoms
->resinfo
[atoms
->atom
[i
].resind
].ic
) &&
75 (std::strcmp(*atoms
->resinfo
[atoms
->atom
[i
].resind
].name
, rresnm
) == 0) &&
76 (std::strstr(*atoms
->atomname
[i
], atomnm
) != nullptr))
83 fprintf(stderr
, "\rCan not find %s%d-%s in pdbfile\n",
84 rresnm
, ri
->nr
, atomnm
);
89 return atoms
->pdbinfo
[i
].bfac
;
92 void correlate_aniso(const char *fn
, t_atoms
*ref
, t_atoms
*calc
,
93 const gmx_output_env_t
*oenv
)
98 fp
= xvgropen(fn
, "Correlation between X-Ray and Computed Uij", "X-Ray",
100 for (i
= 0; (i
< ref
->nr
); i
++)
102 if (ref
->pdbinfo
[i
].bAnisotropic
)
104 for (j
= U11
; (j
<= U23
); j
++)
106 fprintf(fp
, "%10d %10d\n", ref
->pdbinfo
[i
].uij
[j
], calc
->pdbinfo
[i
].uij
[j
]);
113 static void average_residues(double f
[], double **U
, int uind
,
114 int isize
, int index
[], real w_rls
[],
115 const t_atoms
*atoms
)
127 for (i
= 0; i
< isize
; i
++)
129 av
+= w_rls
[index
[i
]]*(f
!= nullptr ? f
[i
] : U
[i
][uind
]);
130 m
+= w_rls
[index
[i
]];
132 atoms
->atom
[index
[i
]].resind
!= atoms
->atom
[index
[i
+1]].resind
)
137 for (j
= start
; j
<= i
; j
++)
144 for (j
= start
; j
<= i
; j
++)
156 void print_dir(FILE *fp
, real
*Uaver
)
158 real eigvec
[DIM
*DIM
];
163 fprintf(fp
, "MSF X Y Z\n");
164 for (d
= 0; d
< DIM
; d
++)
166 fprintf(fp
, " %c ", 'X'+d
-XX
);
167 for (m
= 0; m
< DIM
; m
++)
169 fprintf(fp
, " %9.2e", Uaver
[3*m
+d
]);
171 fprintf(fp
, "%s\n", m
== DIM
? " (nm^2)" : "");
174 for (m
= 0; m
< DIM
*DIM
; m
++)
180 eigensolver(tmp
, DIM
, 0, DIM
, eigval
, eigvec
);
182 fprintf(fp
, "\n Eigenvectors\n\n");
183 fprintf(fp
, "Eigv %-8.2e %-8.2e %-8.2e (nm^2)\n\n",
184 eigval
[2], eigval
[1], eigval
[0]);
185 for (d
= 0; d
< DIM
; d
++)
187 fprintf(fp
, " %c ", 'X'+d
-XX
);
188 for (m
= DIM
-1; m
>= 0; m
--)
190 fprintf(fp
, "%7.4f ", eigvec
[3*m
+d
]);
196 int gmx_rmsf(int argc
, char *argv
[])
198 const char *desc
[] = {
199 "[THISMODULE] computes the root mean square fluctuation (RMSF, i.e. standard ",
200 "deviation) of atomic positions in the trajectory (supplied with [TT]-f[tt])",
201 "after (optionally) fitting to a reference frame (supplied with [TT]-s[tt]).[PAR]",
202 "With option [TT]-oq[tt] the RMSF values are converted to B-factor",
203 "values, which are written to a [REF].pdb[ref] file with the coordinates, of the",
204 "structure file, or of a [REF].pdb[ref] file when [TT]-q[tt] is specified.",
205 "Option [TT]-ox[tt] writes the B-factors to a file with the average",
207 "With the option [TT]-od[tt] the root mean square deviation with",
208 "respect to the reference structure is calculated.[PAR]",
209 "With the option [TT]-aniso[tt], [THISMODULE] will compute anisotropic",
210 "temperature factors and then it will also output average coordinates",
211 "and a [REF].pdb[ref] file with ANISOU records (corresonding to the [TT]-oq[tt]",
212 "or [TT]-ox[tt] option). Please note that the U values",
213 "are orientation-dependent, so before comparison with experimental data",
214 "you should verify that you fit to the experimental coordinates.[PAR]",
215 "When a [REF].pdb[ref] input file is passed to the program and the [TT]-aniso[tt]",
217 "a correlation plot of the Uij will be created, if any anisotropic",
218 "temperature factors are present in the [REF].pdb[ref] file.[PAR]",
219 "With option [TT]-dir[tt] the average MSF (3x3) matrix is diagonalized.",
220 "This shows the directions in which the atoms fluctuate the most and",
223 static gmx_bool bRes
= FALSE
, bAniso
= FALSE
, bFit
= TRUE
;
225 { "-res", FALSE
, etBOOL
, {&bRes
},
226 "Calculate averages for each residue" },
227 { "-aniso", FALSE
, etBOOL
, {&bAniso
},
228 "Compute anisotropic termperature factors" },
229 { "-fit", FALSE
, etBOOL
, {&bFit
},
230 "Do a least squares superposition before computing RMSF. Without this you must make sure that the reference structure and the trajectory match." }
233 int i
, m
, teller
= 0;
238 t_atoms
*pdbatoms
, *refatoms
;
241 rvec
*x
, *pdbx
, *xref
;
245 FILE *fp
; /* the graphics file */
246 const char *devfn
, *dirfn
;
254 real bfac
, pdb_bfac
, *Uaver
;
257 rvec
*rmsd_x
= nullptr;
258 double *rmsf
, invcount
, totmass
;
262 gmx_rmpbc_t gpbc
= nullptr;
264 gmx_output_env_t
*oenv
;
266 const char *leg
[2] = { "MD", "X-Ray" };
269 { efTRX
, "-f", nullptr, ffREAD
},
270 { efTPS
, nullptr, nullptr, ffREAD
},
271 { efNDX
, nullptr, nullptr, ffOPTRD
},
272 { efPDB
, "-q", nullptr, ffOPTRD
},
273 { efPDB
, "-oq", "bfac", ffOPTWR
},
274 { efPDB
, "-ox", "xaver", ffOPTWR
},
275 { efXVG
, "-o", "rmsf", ffWRITE
},
276 { efXVG
, "-od", "rmsdev", ffOPTWR
},
277 { efXVG
, "-oc", "correl", ffOPTWR
},
278 { efLOG
, "-dir", "rmsf", ffOPTWR
}
280 #define NFILE asize(fnm)
282 if (!parse_common_args(&argc
, argv
, PCA_CAN_TIME
| PCA_CAN_VIEW
,
283 NFILE
, fnm
, asize(pargs
), pargs
, asize(desc
), desc
, 0, nullptr,
289 bReadPDB
= ftp2bSet(efPDB
, NFILE
, fnm
);
290 devfn
= opt2fn_null("-od", NFILE
, fnm
);
291 dirfn
= opt2fn_null("-dir", NFILE
, fnm
);
293 read_tps_conf(ftp2fn(efTPS
, NFILE
, fnm
), &top
, &ePBC
, &xref
, nullptr, box
, TRUE
);
294 const char *title
= *top
.name
;
295 snew(w_rls
, top
.atoms
.nr
);
297 fprintf(stderr
, "Select group(s) for root mean square calculation\n");
298 get_index(&top
.atoms
, ftp2fn_null(efNDX
, NFILE
, fnm
), 1, &isize
, &index
, &grpnames
);
301 for (i
= 0; i
< isize
; i
++)
303 w_rls
[index
[i
]] = top
.atoms
.atom
[index
[i
]].m
;
306 /* Malloc the rmsf arrays */
307 snew(xav
, isize
*DIM
);
309 for (i
= 0; i
< isize
; i
++)
323 /* Read coordinates twice */
324 read_tps_conf(opt2fn("-q", NFILE
, fnm
), top_pdb
, nullptr, nullptr, nullptr, pdbbox
, FALSE
);
326 *pdbatoms
= top_pdb
->atoms
;
327 read_tps_conf(opt2fn("-q", NFILE
, fnm
), top_pdb
, nullptr, &pdbx
, nullptr, pdbbox
, FALSE
);
328 /* TODO Should this assert that top_pdb->atoms.nr == top.atoms.nr?
329 * See discussion at https://gerrit.gromacs.org/#/c/6430/1 */
330 title
= *top_pdb
->name
;
332 *refatoms
= top_pdb
->atoms
;
337 pdbatoms
= &top
.atoms
;
338 refatoms
= &top
.atoms
;
340 snew(pdbatoms
->pdbinfo
, pdbatoms
->nr
);
341 pdbatoms
->havePdbInfo
= TRUE
;
342 copy_mat(box
, pdbbox
);
347 sub_xcm(xref
, isize
, index
, top
.atoms
.atom
, xcm
, FALSE
);
350 natom
= read_first_x(oenv
, &status
, ftp2fn(efTRX
, NFILE
, fnm
), &t
, &x
, box
);
354 gpbc
= gmx_rmpbc_init(&top
.idef
, ePBC
, natom
);
357 /* Now read the trj again to compute fluctuations */
363 /* Remove periodic boundary */
364 gmx_rmpbc(gpbc
, natom
, box
, x
);
366 /* Set center of mass to zero */
367 sub_xcm(x
, isize
, index
, top
.atoms
.atom
, xcm
, FALSE
);
369 /* Fit to reference structure */
370 do_fit(natom
, w_rls
, xref
, x
);
373 /* Calculate Anisotropic U Tensor */
374 for (i
= 0; i
< isize
; i
++)
377 for (d
= 0; d
< DIM
; d
++)
379 xav
[i
*DIM
+ d
] += x
[aid
][d
];
380 for (m
= 0; m
< DIM
; m
++)
382 U
[i
][d
*DIM
+ m
] += x
[aid
][d
]*x
[aid
][m
];
389 /* Calculate RMS Deviation */
390 for (i
= 0; (i
< isize
); i
++)
393 for (d
= 0; (d
< DIM
); d
++)
395 rmsd_x
[i
][d
] += gmx::square(x
[aid
][d
]-xref
[aid
][d
]);
402 while (read_next_x(oenv
, status
, &t
, x
, box
));
407 gmx_rmpbc_done(gpbc
);
411 invcount
= 1.0/count
;
412 snew(Uaver
, DIM
*DIM
);
414 for (i
= 0; i
< isize
; i
++)
416 for (d
= 0; d
< DIM
; d
++)
418 xav
[i
*DIM
+ d
] *= invcount
;
420 for (d
= 0; d
< DIM
; d
++)
422 for (m
= 0; m
< DIM
; m
++)
424 U
[i
][d
*DIM
+ m
] = U
[i
][d
*DIM
+ m
]*invcount
425 - xav
[i
*DIM
+ d
]*xav
[i
*DIM
+ m
];
426 Uaver
[3*d
+m
] += top
.atoms
.atom
[index
[i
]].m
*U
[i
][d
*DIM
+ m
];
429 totmass
+= top
.atoms
.atom
[index
[i
]].m
;
431 for (d
= 0; d
< DIM
*DIM
; d
++)
438 for (d
= 0; d
< DIM
*DIM
; d
++)
440 average_residues(nullptr, U
, d
, isize
, index
, w_rls
, &top
.atoms
);
446 for (i
= 0; i
< isize
; i
++)
449 pdbatoms
->pdbinfo
[aid
].bAnisotropic
= TRUE
;
450 pdbatoms
->pdbinfo
[aid
].uij
[U11
] = static_cast<int>(1e6
*U
[i
][XX
*DIM
+ XX
]);
451 pdbatoms
->pdbinfo
[aid
].uij
[U22
] = static_cast<int>(1e6
*U
[i
][YY
*DIM
+ YY
]);
452 pdbatoms
->pdbinfo
[aid
].uij
[U33
] = static_cast<int>(1e6
*U
[i
][ZZ
*DIM
+ ZZ
]);
453 pdbatoms
->pdbinfo
[aid
].uij
[U12
] = static_cast<int>(1e6
*U
[i
][XX
*DIM
+ YY
]);
454 pdbatoms
->pdbinfo
[aid
].uij
[U13
] = static_cast<int>(1e6
*U
[i
][XX
*DIM
+ ZZ
]);
455 pdbatoms
->pdbinfo
[aid
].uij
[U23
] = static_cast<int>(1e6
*U
[i
][YY
*DIM
+ ZZ
]);
467 for (i
= 0; i
< isize
; i
++)
469 rmsf
[i
] = U
[i
][XX
*DIM
+ XX
] + U
[i
][YY
*DIM
+ YY
] + U
[i
][ZZ
*DIM
+ ZZ
];
474 fprintf(stdout
, "\n");
475 print_dir(stdout
, Uaver
);
476 fp
= gmx_ffopen(dirfn
, "w");
477 print_dir(fp
, Uaver
);
481 for (i
= 0; i
< isize
; i
++)
487 /* Write RMSF output */
490 bfac
= 8.0*M_PI
*M_PI
/3.0*100;
491 fp
= xvgropen(ftp2fn(efXVG
, NFILE
, fnm
), "B-Factors",
492 label
, "(A\\b\\S\\So\\N\\S2\\N)", oenv
);
493 xvgr_legend(fp
, 2, leg
, oenv
);
494 for (i
= 0; (i
< isize
); i
++)
496 if (!bRes
|| i
+1 == isize
||
497 top
.atoms
.atom
[index
[i
]].resind
!= top
.atoms
.atom
[index
[i
+1]].resind
)
499 resind
= top
.atoms
.atom
[index
[i
]].resind
;
500 pdb_bfac
= find_pdb_bfac(pdbatoms
, &top
.atoms
.resinfo
[resind
],
501 *(top
.atoms
.atomname
[index
[i
]]));
503 fprintf(fp
, "%5d %10.5f %10.5f\n",
504 bRes
? top
.atoms
.resinfo
[top
.atoms
.atom
[index
[i
]].resind
].nr
: index
[i
]+1, rmsf
[i
]*bfac
,
512 fp
= xvgropen(ftp2fn(efXVG
, NFILE
, fnm
), "RMS fluctuation", label
, "(nm)", oenv
);
513 for (i
= 0; i
< isize
; i
++)
515 if (!bRes
|| i
+1 == isize
||
516 top
.atoms
.atom
[index
[i
]].resind
!= top
.atoms
.atom
[index
[i
+1]].resind
)
518 fprintf(fp
, "%5d %8.4f\n",
519 bRes
? top
.atoms
.resinfo
[top
.atoms
.atom
[index
[i
]].resind
].nr
: index
[i
]+1, std::sqrt(rmsf
[i
]));
525 for (i
= 0; i
< isize
; i
++)
527 pdbatoms
->pdbinfo
[index
[i
]].bfac
= 800*M_PI
*M_PI
/3.0*rmsf
[i
];
532 for (i
= 0; i
< isize
; i
++)
534 rmsf
[i
] = (rmsd_x
[i
][XX
]+rmsd_x
[i
][YY
]+rmsd_x
[i
][ZZ
])/count
;
538 average_residues(rmsf
, nullptr, 0, isize
, index
, w_rls
, &top
.atoms
);
540 /* Write RMSD output */
541 fp
= xvgropen(devfn
, "RMS Deviation", label
, "(nm)", oenv
);
542 for (i
= 0; i
< isize
; i
++)
544 if (!bRes
|| i
+1 == isize
||
545 top
.atoms
.atom
[index
[i
]].resind
!= top
.atoms
.atom
[index
[i
+1]].resind
)
547 fprintf(fp
, "%5d %8.4f\n",
548 bRes
? top
.atoms
.resinfo
[top
.atoms
.atom
[index
[i
]].resind
].nr
: index
[i
]+1, std::sqrt(rmsf
[i
]));
554 if (opt2bSet("-oq", NFILE
, fnm
))
556 /* Write a .pdb file with B-factors and optionally anisou records */
557 for (i
= 0; i
< isize
; i
++)
559 rvec_inc(pdbx
[index
[i
]], xcm
);
561 write_sto_conf_indexed(opt2fn("-oq", NFILE
, fnm
), title
, pdbatoms
, pdbx
,
562 nullptr, ePBC
, pdbbox
, isize
, index
);
564 if (opt2bSet("-ox", NFILE
, fnm
))
567 snew(bFactorX
, top
.atoms
.nr
);
568 for (i
= 0; i
< isize
; i
++)
570 for (d
= 0; d
< DIM
; d
++)
572 bFactorX
[index
[i
]][d
] = xcm
[d
] + xav
[i
*DIM
+ d
];
575 /* Write a .pdb file with B-factors and optionally anisou records */
576 write_sto_conf_indexed(opt2fn("-ox", NFILE
, fnm
), title
, pdbatoms
, bFactorX
, nullptr,
577 ePBC
, pdbbox
, isize
, index
);
582 correlate_aniso(opt2fn("-oc", NFILE
, fnm
), refatoms
, pdbatoms
, oenv
);
583 do_view(oenv
, opt2fn("-oc", NFILE
, fnm
), "-nxy");
585 do_view(oenv
, opt2fn("-o", NFILE
, fnm
), "-nxy");
588 do_view(oenv
, opt2fn("-od", NFILE
, fnm
), "-nxy");