From b0f4bf1d0a83e6bf095d3e557d3a160fe9fac8fc Mon Sep 17 00:00:00 2001 From: Paul Bauer Date: Tue, 22 May 2018 17:01:20 +0200 Subject: [PATCH] Fix PQR file output PQR files from editconf were always written as fixed format PDB files with just the field information added. As pointed out in the linked redmine, this can violate the PQR file standard if the field lengths are too long, even though the file would still be a valid PDB. This adds a slightly different form of the writeout that has a flexible, PQR conform format. Fixes #2511 Change-Id: I626380b642a0214970753da289e9c969ce411ea7 --- docs/release-notes/2018/2018.2.rst | 9 +++ src/gromacs/fileio/confio.cpp | 4 +- src/gromacs/fileio/pdbio.cpp | 110 ++++++++++++++++++++++++++++-------- src/gromacs/fileio/pdbio.h | 5 +- src/gromacs/fileio/trxio.cpp | 2 +- src/gromacs/gmxana/gmx_do_dssp.cpp | 4 +- src/gromacs/gmxana/gmx_editconf.cpp | 17 +++++- 7 files changed, 116 insertions(+), 35 deletions(-) diff --git a/docs/release-notes/2018/2018.2.rst b/docs/release-notes/2018/2018.2.rst index c168c76bfd..139212485d 100644 --- a/docs/release-notes/2018/2018.2.rst +++ b/docs/release-notes/2018/2018.2.rst @@ -40,6 +40,15 @@ Fixed enemat when the .edr file had no matching energy groups :issue:`2508` +Fix PQR file output +"""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""""" +PQR files from gmx editconf violated the standard for the format because +they were always written in fixed format. This commit fixes the issue by +introducing a different output method for PQR files that follows the +standard. + +:issue:`2511` + Fixes to improve portability ^^^^^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/src/gromacs/fileio/confio.cpp b/src/gromacs/fileio/confio.cpp index b9015383b4..cfb8d0fac3 100644 --- a/src/gromacs/fileio/confio.cpp +++ b/src/gromacs/fileio/confio.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -101,7 +101,7 @@ void write_sto_conf_indexed(const char *outfile, const char *title, case efENT: case efPQR: out = gmx_fio_fopen(outfile, "w"); - write_pdbfile_indexed(out, title, atoms, x, ePBC, box, ' ', -1, nindex, index, nullptr, TRUE); + write_pdbfile_indexed(out, title, atoms, x, ePBC, box, ' ', -1, nindex, index, nullptr, TRUE, TRUE); gmx_fio_fclose(out); break; case efESP: diff --git a/src/gromacs/fileio/pdbio.cpp b/src/gromacs/fileio/pdbio.cpp index d6c7579199..da9893d6b4 100644 --- a/src/gromacs/fileio/pdbio.cpp +++ b/src/gromacs/fileio/pdbio.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -259,11 +259,56 @@ static void read_cryst1(char *line, int *ePBC, matrix box) } } +static int +gmx_fprintf_pqr_atomline(FILE * fp, + enum PDB_record record, + int atom_seq_number, + const char * atom_name, + const char * res_name, + char chain_id, + int res_seq_number, + real x, + real y, + real z, + real occupancy, + real b_factor) +{ + GMX_RELEASE_ASSERT(record == epdbATOM || record == epdbHETATM, + "Can only print PQR atom lines as ATOM or HETATM records"); + + /* Check atom name */ + GMX_RELEASE_ASSERT(atom_name != nullptr, + "Need atom information to print pqr"); + + /* Check residue name */ + GMX_RELEASE_ASSERT(res_name != nullptr, + "Need residue information to print pqr"); + + /* Truncate integers so they fit */ + atom_seq_number = atom_seq_number % 100000; + res_seq_number = res_seq_number % 10000; + + int n = fprintf(fp, + "%s %d %s %s %c %d %8.3f %8.3f %8.3f %6.2f %6.2f\n", + pdbtp[record], + atom_seq_number, + atom_name, + res_name, + chain_id, + res_seq_number, + x, y, z, + occupancy, + b_factor); + + return n; +} + void write_pdbfile_indexed(FILE *out, const char *title, const t_atoms *atoms, const rvec x[], int ePBC, const matrix box, char chainid, int model_nr, int nindex, const int index[], - gmx_conect conect, gmx_bool bTerSepChains) + gmx_conect conect, gmx_bool bTerSepChains, + bool usePqrFormat) { gmx_conect_t *gc = (gmx_conect_t *)conect; char resnm[6], nm[6]; @@ -370,29 +415,44 @@ void write_pdbfile_indexed(FILE *out, const char *title, } occup = bOccup ? 1.0 : pdbinfo.occup; bfac = pdbinfo.bfac; - - gmx_fprintf_pdb_atomline(out, - type, - i+1, - nm, - altloc, - resnm, - ch, - resnr, - resic, - 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], - occup, - bfac, - atoms->atom[i].elem); - - if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic) + if (!usePqrFormat) + { + gmx_fprintf_pdb_atomline(out, + type, + i+1, + nm, + altloc, + resnm, + ch, + resnr, + resic, + 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], + occup, + bfac, + atoms->atom[i].elem); + + if (atoms->pdbinfo && atoms->pdbinfo[i].bAnisotropic) + { + fprintf(out, "ANISOU%5d %-4.4s%4.4s%c%4d%c %7d%7d%7d%7d%7d%7d\n", + (i+1)%100000, nm, resnm, ch, resnr, + (resic == '\0') ? ' ' : resic, + atoms->pdbinfo[i].uij[0], atoms->pdbinfo[i].uij[1], + atoms->pdbinfo[i].uij[2], atoms->pdbinfo[i].uij[3], + atoms->pdbinfo[i].uij[4], atoms->pdbinfo[i].uij[5]); + } + } + else { - fprintf(out, "ANISOU%5d %-4.4s%4.4s%c%4d%c %7d%7d%7d%7d%7d%7d\n", - (i+1)%100000, nm, resnm, ch, resnr, - (resic == '\0') ? ' ' : resic, - atoms->pdbinfo[i].uij[0], atoms->pdbinfo[i].uij[1], - atoms->pdbinfo[i].uij[2], atoms->pdbinfo[i].uij[3], - atoms->pdbinfo[i].uij[4], atoms->pdbinfo[i].uij[5]); + gmx_fprintf_pqr_atomline(out, + type, + i+1, + nm, + resnm, + ch, + resnr, + 10*x[i][XX], 10*x[i][YY], 10*x[i][ZZ], + occup, + bfac); } } @@ -422,7 +482,7 @@ void write_pdbfile(FILE *out, const char *title, const t_atoms *atoms, const rve index[i] = i; } write_pdbfile_indexed(out, title, atoms, x, ePBC, box, chainid, model_nr, - atoms->nr, index, conect, bTerSepChains); + atoms->nr, index, conect, bTerSepChains, false); sfree(index); } diff --git a/src/gromacs/fileio/pdbio.h b/src/gromacs/fileio/pdbio.h index 27d5106d53..c88a5c2834 100644 --- a/src/gromacs/fileio/pdbio.h +++ b/src/gromacs/fileio/pdbio.h @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -95,7 +95,8 @@ void gmx_write_pdb_box(FILE *out, int ePBC, const matrix box); void write_pdbfile_indexed(FILE *out, const char *title, const t_atoms *atoms, const rvec x[], int ePBC, const matrix box, char chain, int model_nr, int nindex, const int index[], - gmx_conect conect, gmx_bool bTerSepChains); + gmx_conect conect, gmx_bool bTerSepChains, + bool usePqrFormat); /* REALLY low level */ void write_pdbfile(FILE *out, const char *title, const t_atoms *atoms, diff --git a/src/gromacs/fileio/trxio.cpp b/src/gromacs/fileio/trxio.cpp index a9932699c4..2ed5c41cc1 100644 --- a/src/gromacs/fileio/trxio.cpp +++ b/src/gromacs/fileio/trxio.cpp @@ -427,7 +427,7 @@ int write_trxframe_indexed(t_trxstatus *status, const t_trxframe *fr, int nind, else { write_pdbfile_indexed(gmx_fio_getfp(status->fio), title, fr->atoms, - fr->x, -1, fr->box, ' ', fr->step, nind, ind, gc, TRUE); + fr->x, -1, fr->box, ' ', fr->step, nind, ind, gc, TRUE, FALSE); } break; case efG96: diff --git a/src/gromacs/gmxana/gmx_do_dssp.cpp b/src/gromacs/gmxana/gmx_do_dssp.cpp index a8aab1bc21..039e963881 100644 --- a/src/gromacs/gmxana/gmx_do_dssp.cpp +++ b/src/gromacs/gmxana/gmx_do_dssp.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2012,2013,2014,2015,2017, by the GROMACS development team, led by + * Copyright (c) 2012,2013,2014,2015,2017,2018, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -660,7 +660,7 @@ int gmx_do_dssp(int argc, char *argv[]) } gmx_rmpbc(gpbc, natoms, box, x); tapein = gmx_ffopen(pdbfile, "w"); - write_pdbfile_indexed(tapein, nullptr, atoms, x, ePBC, box, ' ', -1, gnx, index, nullptr, TRUE); + write_pdbfile_indexed(tapein, nullptr, atoms, x, ePBC, box, ' ', -1, gnx, index, nullptr, TRUE, FALSE); gmx_ffclose(tapein); if (0 != system(dssp)) diff --git a/src/gromacs/gmxana/gmx_editconf.cpp b/src/gromacs/gmxana/gmx_editconf.cpp index 7abb509a26..f26bd749cb 100644 --- a/src/gromacs/gmxana/gmx_editconf.cpp +++ b/src/gromacs/gmxana/gmx_editconf.cpp @@ -3,7 +3,7 @@ * * Copyright (c) 1991-2000, University of Groningen, The Netherlands. * Copyright (c) 2001-2004, The GROMACS development team. - * Copyright (c) 2013,2014,2015,2016,2017, by the GROMACS development team, led by + * Copyright (c) 2013,2014,2015,2016,2017,2018, by the GROMACS development team, led by * Mark Abraham, David van der Spoel, Berk Hess, and Erik Lindahl, * and including many others, as listed in the AUTHORS file in the * top-level source directory and at http://www.gromacs.org. @@ -1265,7 +1265,7 @@ int gmx_editconf(int argc, char *argv[]) if (outftp == efPDB) { out = gmx_ffopen(outfile, "w"); - write_pdbfile_indexed(out, *top_tmp->name, &atoms, x, ePBC, box, ' ', 1, isize, index, conect, TRUE); + write_pdbfile_indexed(out, *top_tmp->name, &atoms, x, ePBC, box, ' ', 1, isize, index, conect, TRUE, FALSE); gmx_ffclose(out); } else @@ -1311,7 +1311,18 @@ int gmx_editconf(int argc, char *argv[]) atoms.resinfo[atoms.atom[i].resind].chainid = label[0]; } } - write_pdbfile(out, *top_tmp->name, &atoms, x, ePBC, box, ' ', -1, conect, TRUE); + /* Need to bypass the regular write_pdbfile because I don't want to change + * all instances to include the boolean flag for writing out PQR files. + */ + int *index; + snew(index, atoms.nr); + for (int i = 0; i < atoms.nr; i++) + { + index[i] = i; + } + write_pdbfile_indexed(out, *top_tmp->name, &atoms, x, ePBC, box, ' ', -1, atoms.nr, index, conect, + TRUE, outftp == efPQR ? true : false); + sfree(index); if (bLegend) { pdb_legend(out, atoms.nr, atoms.nres, &atoms, x); -- 2.11.4.GIT