int
gmx_bar(int argc, char *argv[]);
-int
-gmx_bond(int argc, char *argv[]);
-
int
gmx_bundle(int argc, char *argv[]);
int
gmx_disre(int argc, char *argv[]);
-int
-gmx_dist(int argc, char *argv[]);
-
int
gmx_do_dssp(int argc, char *argv[]);
int
gmx_sas(int argc, char *argv[]);
-int
-gmx_select(int argc, char *argv[]);
-
-int
-gmx_sgangle(int argc, char *argv[]);
-
int
gmx_sham(int argc, char *argv[]);
+++ /dev/null
-/*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
- *
- * And Hey:
- * Green Red Orange Magenta Azure Cyan Skyblue
- */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-#include <math.h>
-#include <string.h>
-#include "sysstuff.h"
-#include "typedefs.h"
-#include "smalloc.h"
-#include "macros.h"
-#include "vec.h"
-#include "pbc.h"
-#include "xvgr.h"
-#include "gmx_fatal.h"
-#include "futil.h"
-#include "statutil.h"
-#include "index.h"
-#include "gmx_statistics.h"
-#include "tpxio.h"
-#include "gmx_ana.h"
-
-#ifndef even
-#define even(a) ( ( (a+1) / 2) == (a / 2) )
-#endif
-
-static void make_dist_leg(FILE *fp, int gnx, atom_id index[], t_atoms *atoms,
- const output_env_t oenv)
-{
- char **leg;
- int i;
-
- snew(leg, gnx/2);
- for (i = 0; i < gnx; i += 2)
- {
- snew(leg[i/2], 256);
- sprintf(leg[i/2], "%s %s%d - %s %s%d",
- *(atoms->atomname[index[i]]),
- *(atoms->resinfo[atoms->atom[index[i]].resind].name),
- atoms->resinfo[atoms->atom[index[i]].resind].nr,
- *(atoms->atomname[index[i+1]]),
- *(atoms->resinfo[atoms->atom[index[i+1]].resind].name),
- atoms->resinfo[atoms->atom[index[i+1]].resind].nr);
- }
- xvgr_legend(fp, gnx/2, (const char**)leg, oenv);
- for (i = 0; i < gnx/2; i++)
- {
- sfree(leg[i]);
- }
- sfree(leg);
-}
-
-static void do_bonds(FILE *log, const char *fn, const char *fbond,
- const char *fdist, int gnx, atom_id index[],
- real blen, real tol, gmx_bool bAver,
- t_topology *top, int ePBC, gmx_bool bAverDist,
- const output_env_t oenv)
-{
-#define MAXTAB 1000
- FILE *out, *outd = NULL;
- int *btab = NULL;
- real b0 = 0, b1, db = 0;
- real bond, bav;
- gmx_stats_t b_one = NULL, *b_all = NULL;
- /*real mean, mean2, sqrdev2, sigma2;
- int counter;*/
- rvec *x;
- rvec dx;
- t_trxstatus *status;
- int natoms;
- matrix box;
- real t, fac;
- int bind, i, nframes, i0, i1;
- t_pbc pbc;
- int N;
- real aver, sigma, error;
-
- if (!bAver)
- {
- snew(b_all, gnx/2);
- for (i = 0; (i < gnx/2); i++)
- {
- b_all[i] = gmx_stats_init();
- }
- }
- else
- {
- b_one = gmx_stats_init();
- snew(btab, MAXTAB+1);
- }
-
- natoms = read_first_x(oenv, &status, fn, &t, &x, box);
- if (natoms == 0)
- {
- gmx_fatal(FARGS, "No atoms in trajectory!");
- }
-
- if (fdist)
- {
- outd = xvgropen(fdist, bAverDist ? "Average distance" : "Distances",
- "Time (ps)", "Distance (nm)", oenv);
- if (!bAverDist)
- {
- make_dist_leg(outd, gnx, index, &(top->atoms), oenv);
- }
- }
-
- nframes = 0;
- do
- {
- set_pbc(&pbc, ePBC, box);
- if (fdist)
- {
- fprintf(outd, " %8.4f", t);
- }
- nframes++; /* count frames */
- bav = 0.0;
- for (i = 0; (i < gnx); i += 2)
- {
- pbc_dx(&pbc, x[index[i]], x[index[i+1]], dx);
- bond = norm(dx);
- if (bAverDist)
- {
- bav += bond;
- }
- else if (fdist)
- {
- fprintf(outd, " %.3f", bond);
- }
- if (bAver)
- {
- gmx_stats_add_point(b_one, t, bond, 0, 0);
- if (db == 0)
- {
- if (blen == -1)
- {
- b0 = 0;
- b1 = 0.2;
- db = (b1-b0)/MAXTAB;
- }
- else
- {
- b0 = (1.0-tol)*blen;
- b1 = (1.0+tol)*blen;
- db = (2.0*(b1-b0))/MAXTAB;
- }
- }
- bind = (int)((bond-b0)/db+0.5);
- if ((bind >= 0) && (bind <= MAXTAB))
- {
- btab[bind]++;
- }
- else
- {
- /*
- printf("bond: %4d-%4d bond=%10.5e, dx=(%10.5e,%10.5e,%10.5e)\n",
- index[i],index[i+1],bond,dx[XX],dx[YY],dx[ZZ]);
- */
- }
- }
- else
- {
- gmx_stats_add_point(b_all[i/2], t, bond, 0, 0);
- }
- }
- if (bAverDist)
- {
- fprintf(outd, " %.5f", bav*2.0/gnx);
- }
- if (fdist)
- {
- fprintf(outd, "\n");
- }
- }
- while (read_next_x(oenv, status, &t, x, box));
- close_trj(status);
-
- if (fdist)
- {
- ffclose(outd);
- }
-
- /*
- mean = mean / counter;
- mean2 = mean2 / counter;
- sqrdev2 = (mean2 - mean*mean);
- sigma2 = sqrdev2*counter / (counter - 1);
- */
- /* For definitions see "Weet wat je meet" */
- if (bAver)
- {
- printf("\n");
- gmx_stats_get_npoints(b_one, &N);
- printf("Total number of samples : %d\n", N);
- gmx_stats_get_ase(b_one, &aver, &sigma, &error);
- printf("Mean : %g\n", aver);
- printf("Standard deviation of the distribution: %g\n", sigma);
- printf("Standard deviation of the mean : %g\n", error);
- gmx_stats_done(b_one);
- sfree(b_one);
-
- out = xvgropen(fbond, "Bond Stretching Distribution",
- "Bond Length (nm)", "", oenv);
-
- for (i0 = 0; ((i0 < MAXTAB) && (btab[i0] == 0)); i0++)
- {
- ;
- }
- i0 = max(0, i0-1);
- for (i1 = MAXTAB; ((i1 > 0) && (btab[i1] == 0)); i1--)
- {
- ;
- }
- i1 = min(MAXTAB, i1+1);
-
- if (i0 >= i1)
- {
- gmx_fatal(FARGS, "No distribution... (i0 = %d, i1 = %d)? ? ! ! ? !", i0, i1);
- }
-
- fac = 2.0/(nframes*gnx*db);
- for (i = i0; (i <= i1); i++)
- {
- fprintf(out, "%8.5f %8.5f\n", b0+i*db, btab[i]*fac);
- }
- ffclose(out);
- }
- else
- {
- fprintf(log, "%5s %5s %8s %8s\n", "i", "j", "b_aver", "sigma");
- for (i = 0; (i < gnx/2); i++)
- {
- gmx_stats_get_ase(b_all[i], &aver, &sigma, NULL);
- fprintf(log, "%5u %5u %8.5f %8.5f\n", 1+index[2*i], 1+index[2*i+1],
- aver, sigma);
- gmx_stats_done(b_all[i]);
- sfree(b_all[i]);
- }
- sfree(b_all);
- }
-}
-
-int gmx_bond(int argc, char *argv[])
-{
- const char *desc[] = {
- "[TT]g_bond[tt] makes a distribution of bond lengths. If all is well a",
- "Gaussian distribution should be made when using a harmonic potential.",
- "Bonds are read from a single group in the index file in order i1-j1",
- "i2-j2 through in-jn.[PAR]",
- "[TT]-tol[tt] gives the half-width of the distribution as a fraction",
- "of the bondlength ([TT]-blen[tt]). That means, for a bond of 0.2",
- "a tol of 0.1 gives a distribution from 0.18 to 0.22.[PAR]",
- "Option [TT]-d[tt] plots all the distances as a function of time.",
- "This requires a structure file for the atom and residue names in",
- "the output. If however the option [TT]-averdist[tt] is given (as well",
- "or separately) the average bond length is plotted instead."
- };
- const char *bugs[] = {
- "It should be possible to get bond information from the topology."
- };
- static real blen = -1.0, tol = 0.1;
- static gmx_bool bAver = TRUE, bAverDist = TRUE;
- t_pargs pa[] = {
- { "-blen", FALSE, etREAL, {&blen},
- "Bond length. By default length of first bond" },
- { "-tol", FALSE, etREAL, {&tol},
- "Half width of distribution as fraction of [TT]-blen[tt]" },
- { "-aver", FALSE, etBOOL, {&bAver},
- "Average bond length distributions" },
- { "-averdist", FALSE, etBOOL, {&bAverDist},
- "Average distances (turns on [TT]-d[tt])" }
- };
- FILE *fp;
- char *grpname;
- const char *fdist;
- int gnx;
- atom_id *index;
- char title[STRLEN];
- t_topology top;
- int ePBC = -1;
- rvec *x;
- matrix box;
- output_env_t oenv;
-
- t_filenm fnm[] = {
- { efTRX, "-f", NULL, ffREAD },
- { efNDX, NULL, NULL, ffREAD },
- { efTPS, NULL, NULL, ffOPTRD },
- { efXVG, "-o", "bonds", ffWRITE },
- { efLOG, NULL, "bonds", ffOPTWR },
- { efXVG, "-d", "distance", ffOPTWR }
- };
-#define NFILE asize(fnm)
-
- if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
- NFILE, fnm, asize(pa), pa, asize(desc), desc, asize(bugs), bugs,
- &oenv))
- {
- return 0;
- }
-
- if (bAverDist)
- {
- fdist = opt2fn("-d", NFILE, fnm);
- }
- else
- {
- fdist = opt2fn_null("-d", NFILE, fnm);
- if (fdist)
- {
- read_tps_conf(ftp2fn(efTPS, NFILE, fnm), title, &top, &ePBC, &x, NULL, box,
- FALSE);
- }
- }
-
- rd_index(ftp2fn(efNDX, NFILE, fnm), 1, &gnx, &index, &grpname);
- if (!even(gnx) )
- {
- fprintf(stderr, "WARNING: odd number of atoms (%d) in group!\n", gnx);
- }
- fprintf(stderr, "Will gather information on %d bonds\n", gnx/2);
-
- if (!bAver)
- {
- fp = ftp2FILE(efLOG, NFILE, fnm, "w");
- }
- else
- {
- fp = NULL;
- }
-
- do_bonds(fp, ftp2fn(efTRX, NFILE, fnm), opt2fn("-o", NFILE, fnm), fdist, gnx, index,
- blen, tol, bAver, &top, ePBC, bAverDist, oenv);
-
- do_view(oenv, opt2fn("-o", NFILE, fnm), "-nxy");
- do_view(oenv, opt2fn_null("-d", NFILE, fnm), "-nxy");
-
- return 0;
-}
+++ /dev/null
-/*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
- *
- * And Hey:
- * Green Red Orange Magenta Azure Cyan Skyblue
- */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-#include <typedefs.h>
-
-#include "smalloc.h"
-#include "macros.h"
-#include <math.h>
-#include "xvgr.h"
-#include "statutil.h"
-#include "string2.h"
-#include "vec.h"
-#include "index.h"
-#include "pbc.h"
-#include "gmx_fatal.h"
-#include "futil.h"
-#include "gstat.h"
-#include "gmx_ana.h"
-
-
-static void add_contact_time(int **ccount, int *ccount_nalloc, int t)
-{
- int i;
-
- if (t+2 >= *ccount_nalloc)
- {
- srenew(*ccount, t+2);
- for (i = *ccount_nalloc; i < t+2; i++)
- {
- (*ccount)[i] = 0;
- }
- *ccount_nalloc = t+2;
- }
- (*ccount)[t]++;
-}
-
-int gmx_dist(int argc, char *argv[])
-{
- const char *desc[] = {
- "[TT]g_dist[tt] can calculate the distance between the centers of mass of two",
- "groups of atoms as a function of time. The total distance and its",
- "[IT]x[it]-, [IT]y[it]-, and [IT]z[it]-components are plotted.[PAR]",
- "Or when [TT]-dist[tt] is set, print all the atoms in group 2 that are",
- "closer than a certain distance to the center of mass of group 1.[PAR]",
- "With options [TT]-lt[tt] and [TT]-dist[tt] the number of contacts",
- "of all atoms in group 2 that are closer than a certain distance",
- "to the center of mass of group 1 are plotted as a function of the time",
- "that the contact was continuously present. The [TT]-intra[tt] switch enables",
- "calculations of intramolecular distances avoiding distance calculation to its",
- "periodic images. For a proper function, the molecule in the input trajectory",
- "should be whole (e.g. by preprocessing with [TT]trjconv -pbc[tt]) or a matching",
- "topology should be provided. The [TT]-intra[tt] switch will only give",
- "meaningful results for intramolecular and not intermolecular distances.[PAR]",
- "Other programs that calculate distances are [TT]g_mindist[tt]",
- "and [TT]g_bond[tt]."
- };
-
- t_topology *top = NULL;
- int ePBC;
- real t, t0, cut2, dist2;
- rvec *x = NULL, *v = NULL, dx;
- matrix box;
- t_trxstatus *status;
- int natoms;
-
- int g, d, i, j, res, teller = 0;
- atom_id aid;
-
- int ngrps; /* the number of index groups */
- atom_id **index, natoms_groups; /* the index for the atom numbers */
- int *isize; /* the size of each group */
- char **grpname; /* the name of each group */
- rvec *com;
- real *mass;
- FILE *fp = NULL, *fplt = NULL;
- gmx_bool bCutoff, bPrintDist, bLifeTime, bIntra = FALSE;
- t_pbc *pbc;
- int *contact_time = NULL, *ccount = NULL, ccount_nalloc = 0, sum;
- char buf[STRLEN];
- output_env_t oenv;
- gmx_rmpbc_t gpbc = NULL;
-
- const char *leg[4] = { "|d|", "d\\sx\\N", "d\\sy\\N", "d\\sz\\N" };
-
- static real cut = 0;
-
- t_pargs pa[] = {
- { "-intra", FALSE, etBOOL, {&bIntra},
- "Calculate distances without considering periodic boundaries, e.g. intramolecular." },
- { "-dist", FALSE, etREAL, {&cut},
- "Print all atoms in group 2 closer than dist to the center of mass of group 1" }
- };
-#define NPA asize(pa)
-
- t_filenm fnm[] = {
- { efTRX, "-f", NULL, ffREAD },
- { efTPX, NULL, NULL, ffREAD },
- { efNDX, NULL, NULL, ffOPTRD },
- { efXVG, NULL, "dist", ffOPTWR },
- { efXVG, "-lt", "lifetime", ffOPTWR },
- };
-#define NFILE asize(fnm)
-
- if (!parse_common_args(&argc, argv, PCA_CAN_TIME | PCA_BE_NICE,
- NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
- {
- return 0;
- }
-
- bCutoff = opt2parg_bSet("-dist", NPA, pa);
- cut2 = cut*cut;
- bLifeTime = opt2bSet("-lt", NFILE, fnm);
- bPrintDist = (bCutoff && !bLifeTime);
-
- top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC);
-
- /* read index files */
- ngrps = 2;
- snew(com, ngrps);
- snew(grpname, ngrps);
- snew(index, ngrps);
- snew(isize, ngrps);
- get_index(&top->atoms, ftp2fn_null(efNDX, NFILE, fnm), ngrps, isize, index, grpname);
-
- /* calculate mass */
- natoms_groups = 0;
- snew(mass, ngrps);
- for (g = 0; (g < ngrps); g++)
- {
- mass[g] = 0;
- for (i = 0; (i < isize[g]); i++)
- {
- if (index[g][i] > natoms_groups)
- {
- natoms_groups = index[g][i];
- }
- if (index[g][i] >= top->atoms.nr)
- {
- gmx_fatal(FARGS, "Atom number %d, item %d of group %d, is larger than number of atoms in the topolgy (%d)\n", index[g][i]+1, i+1, g+1, top->atoms.nr+1);
- }
- mass[g] += top->atoms.atom[index[g][i]].m;
- }
- }
- /* The number is one more than the highest index */
- natoms_groups++;
-
- natoms = read_first_x(oenv, &status, ftp2fn(efTRX, NFILE, fnm), &t, &x, box);
- t0 = t;
-
- if (natoms_groups > natoms)
- {
- gmx_fatal(FARGS, "Atom number %d in an index group is larger than number of atoms in the trajectory (%d)\n", (int)natoms_groups, natoms);
- }
-
- if (!bCutoff)
- {
- /* open output file */
- fp = xvgropen(ftp2fn(efXVG, NFILE, fnm),
- "Distance", "Time (ps)", "Distance (nm)", oenv);
- xvgr_legend(fp, 4, leg, oenv);
- }
- else
- {
- ngrps = 1;
- if (bLifeTime)
- {
- snew(contact_time, isize[1]);
- }
- }
- if (ePBC != epbcNONE)
- {
- snew(pbc, 1);
- }
- else
- {
- pbc = NULL;
- }
-
- gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
- do
- {
- /* initialisation for correct distance calculations */
- if (pbc)
- {
- set_pbc(pbc, ePBC, box);
- /* make molecules whole again */
- gmx_rmpbc(gpbc, natoms, box, x);
- }
- /* calculate center of masses */
- for (g = 0; (g < ngrps); g++)
- {
- if (isize[g] == 1)
- {
- copy_rvec(x[index[g][0]], com[g]);
- }
- else
- {
- for (d = 0; (d < DIM); d++)
- {
- com[g][d] = 0;
- for (i = 0; (i < isize[g]); i++)
- {
- com[g][d] += x[index[g][i]][d] * top->atoms.atom[index[g][i]].m;
- }
- com[g][d] /= mass[g];
- }
- }
- }
-
- if (!bCutoff)
- {
- /* write to output */
- fprintf(fp, "%12.7f ", t);
- for (g = 0; (g < ngrps/2); g++)
- {
- if (pbc && (!bIntra))
- {
- pbc_dx(pbc, com[2*g], com[2*g+1], dx);
- }
- else
- {
- rvec_sub(com[2*g], com[2*g+1], dx);
- }
-
- fprintf(fp, "%12.7f %12.7f %12.7f %12.7f",
- norm(dx), dx[XX], dx[YY], dx[ZZ]);
- }
- fprintf(fp, "\n");
- }
- else
- {
- for (i = 0; (i < isize[1]); i++)
- {
- j = index[1][i];
- if (pbc && (!bIntra))
- {
- pbc_dx(pbc, x[j], com[0], dx);
- }
- else
- {
- rvec_sub(x[j], com[0], dx);
- }
-
- dist2 = norm2(dx);
- if (dist2 < cut2)
- {
- if (bPrintDist)
- {
- res = top->atoms.atom[j].resind;
- fprintf(stdout, "\rt: %g %d %s %d %s %g (nm)\n",
- t, top->atoms.resinfo[res].nr, *top->atoms.resinfo[res].name,
- j+1, *top->atoms.atomname[j], sqrt(dist2));
- }
- if (bLifeTime)
- {
- contact_time[i]++;
- }
- }
- else
- {
- if (bLifeTime)
- {
- if (contact_time[i])
- {
- add_contact_time(&ccount, &ccount_nalloc, contact_time[i]-1);
- contact_time[i] = 0;
- }
- }
- }
- }
- }
-
- teller++;
- }
- while (read_next_x(oenv, status, &t, x, box));
- gmx_rmpbc_done(gpbc);
-
- if (!bCutoff)
- {
- ffclose(fp);
- }
-
- close_trj(status);
-
- if (bCutoff && bLifeTime)
- {
- /* Add the contacts still present in the last frame */
- for (i = 0; i < isize[1]; i++)
- {
- if (contact_time[i])
- {
- add_contact_time(&ccount, &ccount_nalloc, contact_time[i]-1);
- }
- }
-
- sprintf(buf, "%s - %s within %g nm",
- grpname[0], grpname[1], cut);
- fp = xvgropen(opt2fn("-lt", NFILE, fnm),
- buf, "Time (ps)", "Number of contacts", oenv);
- for (i = 0; i < min(ccount_nalloc, teller-1); i++)
- {
- /* Account for all subintervals of longer intervals */
- sum = 0;
- for (j = i; j < ccount_nalloc; j++)
- {
- sum += (j-i+1)*ccount[j];
- }
-
- fprintf(fp, "%10.3f %10.3f\n", i*(t-t0)/(teller-1), sum/(double)(teller-i));
- }
- ffclose(fp);
- }
-
- return 0;
-}
+++ /dev/null
-/*
- *
- * This source code is part of
- *
- * G R O M A C S
- *
- * GROningen MAchine for Chemical Simulations
- *
- * VERSION 3.2.0
- * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
- * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
- * Copyright (c) 2001-2004, The GROMACS development team,
- * check out http://www.gromacs.org for more information.
-
- * This program is free software; you can redistribute it and/or
- * modify it under the terms of the GNU General Public License
- * as published by the Free Software Foundation; either version 2
- * of the License, or (at your option) any later version.
- *
- * If you want to redistribute modifications, please consider that
- * scientific software is very special. Version control is crucial -
- * bugs must be traceable. We will be happy to consider code for
- * inclusion in the official distribution, but derived work must not
- * be called official GROMACS. Details are found in the README & COPYING
- * files - if they are missing, get the official version at www.gromacs.org.
- *
- * To help us fund GROMACS development, we humbly ask that you cite
- * the papers on the package - you can find them in the top README file.
- *
- * For more info, check our website at http://www.gromacs.org
- *
- * And Hey:
- * Green Red Orange Magenta Azure Cyan Skyblue
- */
-#ifdef HAVE_CONFIG_H
-#include <config.h>
-#endif
-#include <math.h>
-
-#include "sysstuff.h"
-#include <string.h>
-#include "typedefs.h"
-#include "smalloc.h"
-#include "macros.h"
-#include "rmpbc.h"
-#include "vec.h"
-#include "xvgr.h"
-#include "pbc.h"
-#include "futil.h"
-#include "statutil.h"
-#include "index.h"
-#include "gmx_ana.h"
-
-
-/* this version only works correctly if one of the entries in the index file
- is a plane (three atoms specified) and the other a vector. Distance
- is calculated from the center of the plane to both atoms of the vector */
-
-static void print_types(atom_id index1[], int gnx1, char *group1,
- atom_id index2[], int gnx2, char *group2,
- t_topology *top)
-{
- int i, j;
-
- fprintf(stderr, "\n");
- fprintf(stderr, "Group %s contains the following atoms: \n", group1);
- for (i = 0; i < gnx1; i++)
- {
- fprintf(stderr, "Atomname %d: %s\n", i, *(top->atoms.atomname[index1[i]]));
- }
- fprintf(stderr, "\n");
-
- fprintf(stderr, "Group %s contains the following atoms: \n", group2);
- for (j = 0; j < gnx2; j++)
- {
- fprintf(stderr, "Atomname %d: %s\n", j, *(top->atoms.atomname[index2[j]]));
- }
- fprintf(stderr, "\n");
-
- fprintf(stderr, "Careful: distance only makes sense in some situations.\n\n");
-}
-
-static void calculate_normal(atom_id index[], rvec x[], rvec result, rvec center)
-{
- rvec c1, c2;
- int i;
-
- /* calculate centroid of triangle spanned by the three points */
- for (i = 0; i < 3; i++)
- {
- center[i] = (x[index[0]][i] + x[index[1]][i] + x[index[2]][i])/3;
- }
-
- /* use P1P2 x P1P3 to calculate normal, given three points P1-P3 */
- rvec_sub(x[index[1]], x[index[0]], c1); /* find two vectors */
- rvec_sub(x[index[2]], x[index[0]], c2);
-
- cprod(c1, c2, result); /* take crossproduct between these */
-}
-
-/* calculate the angle and distance between the two groups */
-static void calc_angle(int ePBC, matrix box, rvec x[], atom_id index1[],
- atom_id index2[], int gnx1, int gnx2,
- real *angle, real *distance,
- real *distance1, real *distance2)
-
-/* distance is distance between centers, distance 1 between center of plane
- and atom one of vector, distance 2 same for atom two
- */
-
-{
- rvec
- normal1, normal2, /* normals on planes of interest */
- center1, center2, /* center of triangle of points given to define plane,*/
- /* or center of vector if a vector is given */
- h1, h2, h3, h4, h5; /* temp. vectors */
- t_pbc pbc;
-
- set_pbc(&pbc, ePBC, box);
-
- switch (gnx1)
- {
- case 3: /* group 1 defines plane */
- calculate_normal(index1, x, normal1, center1);
- break;
- case 2: /* group 1 defines vector */
- rvec_sub(x[index1[0]], x[index1[1]], normal1);
- rvec_add(x[index1[0]], x[index1[1]], h1);
- svmul(0.5, h1, center1); /* center is geometric mean */
- break;
- default: /* group 1 does none of the above */
- gmx_fatal(FARGS, "Something wrong with contents of index file. Groups should contain 2 or 3 atoms.\n");
- }
-
- switch (gnx2)
- {
- case 3: /* group 2 defines plane */
- calculate_normal(index2, x, normal2, center2);
- break;
- case 2: /* group 2 defines vector */
- rvec_sub(x[index2[0]], x[index2[1]], normal2);
- rvec_add(x[index2[0]], x[index2[1]], h2);
- svmul(0.5, h2, center2); /* center is geometric mean */
- break;
- case 0:
- normal2[XX] = 0;
- normal2[YY] = 0;
- normal2[ZZ] = 1;
- center2[XX] = 0;
- center2[YY] = 0;
- center2[ZZ] = 0;
- break;
- default: /* group 2 does none of the above */
- gmx_fatal(FARGS, "Something wrong with contents of index file.\n");
- }
-
- *angle = cos_angle(normal1, normal2);
-
- if (box)
- {
- pbc_dx(&pbc, center1, center2, h3);
- }
- else
- {
- rvec_sub(center1, center2, h3);
- }
- *distance = norm(h3);
-
- if (gnx1 == 3 && gnx2 == 2)
- {
- if (box)
- {
- pbc_dx(&pbc, center1, x[index2[0]], h4);
- pbc_dx(&pbc, center1, x[index2[1]], h5);
- }
- else
- {
- rvec_sub(center1, x[index2[0]], h4);
- rvec_sub(center1, x[index2[1]], h5);
- }
- *distance1 = norm(h4);
- *distance2 = norm(h5);
- }
- else if (gnx1 == 2 && gnx2 == 3)
- {
- rvec_sub(center1, x[index1[0]], h4);
- rvec_sub(center1, x[index1[1]], h5);
- *distance1 = norm(h4);
- *distance2 = norm(h5);
- }
- else
- {
- *distance1 = 0; *distance2 = 0;
- }
-}
-
-void sgangle_plot(const char *fn, const char *afile, const char *dfile,
- const char *d1file, const char *d2file,
- atom_id index1[], int gnx1, char *grpn1,
- atom_id index2[], int gnx2, char *grpn2,
- t_topology *top, int ePBC, const output_env_t oenv)
-{
- FILE
- *sg_angle, /* xvgr file with angles */
- *sg_distance = NULL, /* xvgr file with distances */
- *sg_distance1 = NULL, /* xvgr file with distance between plane and atom */
- *sg_distance2 = NULL; /* xvgr file with distance between plane and atom2 */
- real
- t, /* time */
- angle, /* cosine of angle between two groups */
- distance, /* distance between two groups. */
- distance1, /* distance between plane and one of two atoms */
- distance2; /* same for second of two atoms */
- t_trxstatus *status;
- int natoms, teller = 0;
- rvec *x0; /* coordinates, and coordinates corrected for pb */
- matrix box;
- char buf[256]; /* for xvgr title */
- gmx_rmpbc_t gpbc = NULL;
- const char* aleg[2] = { "cos(Angle)", "Angle (degrees)" }; /* legends for sg_angle output file */
-
- if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
- {
- gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
- }
-
- sprintf(buf, "Angle between %s and %s", grpn1, grpn2);
- sg_angle = xvgropen(afile, buf, "Time (ps)", "Angle (degrees)", oenv);
- xvgr_legend(sg_angle, 2, aleg, oenv);
-
- if (dfile)
- {
- sprintf(buf, "Distance between %s and %s", grpn1, grpn2);
- sg_distance = xvgropen(dfile, buf, "Time (ps)", "Distance (nm)", oenv);
- }
-
- if (d1file)
- {
- sprintf(buf, "Distance between plane and first atom of vector");
- sg_distance1 = xvgropen(d1file, buf, "Time (ps)", "Distance (nm)", oenv);
- }
-
- if (d2file)
- {
- sprintf(buf, "Distance between plane and second atom of vector");
- sg_distance2 = xvgropen(d2file, buf, "Time (ps", "Distance (nm)", oenv);
- }
-
- gpbc = gmx_rmpbc_init(&(top->idef), ePBC, natoms);
-
- do
- {
- teller++;
-
- gmx_rmpbc(gpbc, natoms, box, x0);
-
- calc_angle(ePBC, box, x0, index1, index2, gnx1, gnx2, &angle,
- &distance, &distance1, &distance2);
-
- fprintf(sg_angle, "%12g %12g %12g\n", t, angle, acos(angle)*180.0/M_PI);
- if (dfile)
- {
- fprintf(sg_distance, "%12g %12g\n", t, distance);
- }
- if (d1file)
- {
- fprintf(sg_distance1, "%12g %12g\n", t, distance1);
- }
- if (d2file)
- {
- fprintf(sg_distance2, "%12g %12g\n", t, distance1);
- }
-
- }
- while (read_next_x(oenv, status, &t, x0, box));
-
- gmx_rmpbc_done(gpbc);
-
- fprintf(stderr, "\n");
- close_trj(status);
- ffclose(sg_angle);
- if (dfile)
- {
- ffclose(sg_distance);
- }
- if (d1file)
- {
- ffclose(sg_distance1);
- }
- if (d2file)
- {
- ffclose(sg_distance2);
- }
-
- sfree(x0);
-}
-
-static void calc_angle_single(int ePBC,
- matrix box,
- rvec xzero[],
- rvec x[],
- atom_id index1[],
- atom_id index2[],
- int gnx1,
- int gnx2,
- real *angle,
- real *distance,
- real *distance1,
- real *distance2)
-{
- t_pbc pbc;
-
- /* distance is distance between centers, distance 1 between center of plane
- and atom one of vector, distance 2 same for atom two
- */
-
- rvec normal1, normal2; /* normals on planes of interest */
- rvec center1, center2;
- /* center of triangle of pts to define plane,
- * or center of vector if a vector is given
- */
- rvec h1, h2, h3, h4, h5; /* temp. vectors */
-
- if (box)
- {
- set_pbc(&pbc, ePBC, box);
- }
-
- switch (gnx1)
- {
- case 3: /* group 1 defines plane */
- calculate_normal(index1, xzero, normal1, center1);
- break;
- case 2: /* group 1 defines vector */
- rvec_sub(xzero[index1[0]], xzero[index1[1]], normal1);
- rvec_add(xzero[index1[0]], xzero[index1[1]], h1);
- svmul(0.5, h1, center1); /* center is geometric mean */
- break;
- default: /* group 1 does none of the above */
- gmx_fatal(FARGS, "Something wrong with contents of index file.\n");
- }
-
- switch (gnx2)
- {
- case 3: /* group 2 defines plane */
- calculate_normal(index2, x, normal2, center2);
- break;
- case 2: /* group 2 defines vector */
- rvec_sub(x[index2[0]], x[index2[1]], normal2);
- rvec_add(x[index2[0]], x[index2[1]], h2);
- svmul(0.5, h2, center2); /* center is geometric mean */
- break;
- default: /* group 2 does none of the above */
- gmx_fatal(FARGS, "Something wrong with contents of index file.\n");
- }
-
- *angle = cos_angle(normal1, normal2);
-
- if (box)
- {
- pbc_dx(&pbc, center1, center2, h3);
- }
- else
- {
- rvec_sub(center1, center2, h3);
- }
- *distance = norm(h3);
-
- if (gnx1 == 3 && gnx2 == 2)
- {
- if (box)
- {
- pbc_dx(&pbc, center1, x[index2[0]], h4);
- pbc_dx(&pbc, center1, x[index2[1]], h5);
- }
- else
- {
- rvec_sub(center1, x[index2[0]], h4);
- rvec_sub(center1, x[index2[1]], h5);
- }
- *distance1 = norm(h4);
- *distance2 = norm(h5);
- }
- else if (gnx1 == 2 && gnx2 == 3)
- {
- rvec_sub(center1, xzero[index1[0]], h4);
- rvec_sub(center1, xzero[index1[1]], h5);
- *distance1 = norm(h4);
- *distance2 = norm(h5);
- }
- else
- {
- *distance1 = 0; *distance2 = 0;
- }
-}
-
-
-void sgangle_plot_single(const char *fn, const char *afile, const char *dfile,
- const char *d1file, const char *d2file,
- atom_id index1[], int gnx1, char *grpn1,
- atom_id index2[], int gnx2, char *grpn2,
- t_topology *top, int ePBC, const output_env_t oenv)
-{
- FILE
- *sg_angle, /* xvgr file with angles */
- *sg_distance = NULL, /* xvgr file with distances */
- *sg_distance1 = NULL, /* xvgr file with distance between plane and atom */
- *sg_distance2 = NULL; /* xvgr file with distance between plane and atom2 */
- real
- t, /* time */
- angle, /* cosine of angle between two groups */
- distance, /* distance between two groups. */
- distance1, /* distance between plane and one of two atoms */
- distance2; /* same for second of two atoms */
- t_trxstatus *status;
- int natoms, teller = 0;
- int i;
- rvec *x0; /* coordinates, and coordinates corrected for pb */
- rvec *xzero;
- matrix box;
- char buf[256]; /* for xvgr title */
- gmx_rmpbc_t gpbc = NULL;
-
-
- if ((natoms = read_first_x(oenv, &status, fn, &t, &x0, box)) == 0)
- {
- gmx_fatal(FARGS, "Could not read coordinates from statusfile\n");
- }
-
- sprintf(buf, "Angle between %s and %s", grpn1, grpn2);
- sg_angle = xvgropen(afile, buf, "Time (ps)", "Cos(angle) ", oenv);
-
- if (dfile)
- {
- sprintf(buf, "Distance between %s and %s", grpn1, grpn2);
- sg_distance = xvgropen(dfile, buf, "Time (ps)", "Distance (nm)", oenv);
- }
-
- if (d1file)
- {
- sprintf(buf, "Distance between plane and first atom of vector");
- sg_distance1 = xvgropen(d1file, buf, "Time (ps)", "Distance (nm)", oenv);
- }
-
- if (d2file)
- {
- sprintf(buf, "Distance between plane and second atom of vector");
- sg_distance2 = xvgropen(d2file, buf, "Time (ps", "Distance (nm)", oenv);
- }
-
- snew(xzero, natoms);
- gpbc = gmx_rmpbc_init(&top->idef, ePBC, natoms);
-
- do
- {
- teller++;
-
- gmx_rmpbc(gpbc, natoms, box, x0);
- if (teller == 1)
- {
- for (i = 0; i < natoms; i++)
- {
- copy_rvec(x0[i], xzero[i]);
- }
- }
-
-
- calc_angle_single(ePBC, box,
- xzero, x0, index1, index2, gnx1, gnx2, &angle,
- &distance, &distance1, &distance2);
-
- fprintf(sg_angle, "%12g %12g %12g\n", t, angle, acos(angle)*180.0/M_PI);
- if (dfile)
- {
- fprintf(sg_distance, "%12g %12g\n", t, distance);
- }
- if (d1file)
- {
- fprintf(sg_distance1, "%12g %12g\n", t, distance1);
- }
- if (d2file)
- {
- fprintf(sg_distance2, "%12g %12g\n", t, distance1);
- }
-
- }
- while (read_next_x(oenv, status, &t, x0, box));
- gmx_rmpbc_done(gpbc);
-
- fprintf(stderr, "\n");
- close_trj(status);
- ffclose(sg_angle);
- if (dfile)
- {
- ffclose(sg_distance);
- }
- if (d1file)
- {
- ffclose(sg_distance1);
- }
- if (d2file)
- {
- ffclose(sg_distance2);
- }
-
- sfree(x0);
-}
-
-
-
-int gmx_sgangle(int argc, char *argv[])
-{
- const char *desc[] = {
- "Compute the angle and distance between two groups. ",
- "The groups are defined by a number of atoms given in an index file and",
- "may be two or three atoms in size.",
- "If [TT]-one[tt] is set, only one group should be specified in the index",
- "file and the angle between this group at time 0 and t will be computed.",
- "The angles calculated depend on the order in which the atoms are ",
- "given. Giving, for instance, 5 6 will rotate the vector 5-6 with ",
- "180 degrees compared to giving 6 5. [PAR]If three atoms are given, ",
- "the normal on the plane spanned by those three atoms will be",
- "calculated, using the formula P1P2 x P1P3.",
- "The cos of the angle is calculated, using the inproduct of the two",
- "normalized vectors.[PAR]",
- "Here is what some of the file options do:[BR]",
- "[TT]-oa[tt]: Angle between the two groups specified in the index file. If a group contains three atoms the normal to the plane defined by those three atoms will be used. If a group contains two atoms, the vector defined by those two atoms will be used.[BR]",
- "[TT]-od[tt]: Distance between two groups. Distance is taken from the center of one group to the center of the other group.[BR]",
- "[TT]-od1[tt]: If one plane and one vector is given, the distances for each of the atoms from the center of the plane is given separately.[BR]",
- "[TT]-od2[tt]: For two planes this option has no meaning."
- };
-
- output_env_t oenv;
- const char *fna, *fnd, *fnd1, *fnd2;
- char * grpname[2]; /* name of the two groups */
- int gnx[2]; /* size of the two groups */
- t_topology *top; /* topology */
- int ePBC;
- atom_id *index[2];
- static gmx_bool bOne = FALSE, bZ = FALSE;
- t_pargs pa[] = {
- { "-one", FALSE, etBOOL, {&bOne},
- "Only one group compute angle between vector at time zero and time t"},
- { "-z", FALSE, etBOOL, {&bZ},
- "Use the [IT]z[it]-axis as reference" }
- };
-#define NPA asize(pa)
-
- t_filenm fnm[] = { /* files for g_sgangle */
- { efTRX, "-f", NULL, ffREAD }, /* trajectory file */
- { efNDX, NULL, NULL, ffREAD }, /* index file */
- { efTPX, NULL, NULL, ffREAD }, /* topology file */
- { efXVG, "-oa", "sg_angle", ffWRITE }, /* xvgr output file */
- { efXVG, "-od", "sg_dist", ffOPTWR }, /* xvgr output file */
- { efXVG, "-od1", "sg_dist1", ffOPTWR }, /* xvgr output file */
- { efXVG, "-od2", "sg_dist2", ffOPTWR } /* xvgr output file */
- };
-
-#define NFILE asize(fnm)
-
- if (!parse_common_args(&argc, argv, PCA_CAN_VIEW | PCA_CAN_TIME | PCA_BE_NICE,
- NFILE, fnm, NPA, pa, asize(desc), desc, 0, NULL, &oenv))
- {
- return 0;
- }
-
-
- top = read_top(ftp2fn(efTPX, NFILE, fnm), &ePBC); /* read topology file */
-
- fna = opt2fn("-oa", NFILE, fnm);
- fnd = opt2fn_null("-od", NFILE, fnm);
- fnd1 = opt2fn_null("-od1", NFILE, fnm);
- fnd2 = opt2fn_null("-od2", NFILE, fnm);
-
- /* read index file. */
- if (bOne)
- {
- rd_index(ftp2fn(efNDX, NFILE, fnm), 1, gnx, index, grpname);
- print_types(index[0], gnx[0], grpname[0],
- index[0], gnx[0], grpname[0], top);
-
- sgangle_plot_single(ftp2fn(efTRX, NFILE, fnm), fna, fnd, fnd1, fnd2,
- index[0], gnx[0], grpname[0],
- index[0], gnx[0], grpname[0],
- top, ePBC, oenv);
- }
- else
- {
- rd_index(ftp2fn(efNDX, NFILE, fnm), bZ ? 1 : 2, gnx, index, grpname);
- if (!bZ)
- {
- print_types(index[0], gnx[0], grpname[0],
- index[1], gnx[1], grpname[1], top);
- }
- else
- {
- gnx[1] = 0;
- grpname[1] = strdup("Z-axis");
- }
- sgangle_plot(ftp2fn(efTRX, NFILE, fnm), fna, fnd, fnd1, fnd2,
- index[0], gnx[0], grpname[0],
- index[1], gnx[1], grpname[1],
- top, ePBC, oenv);
- }
-
- do_view(oenv, fna, "-nxy"); /* view xvgr file */
- do_view(oenv, fnd, "-nxy"); /* view xvgr file */
- do_view(oenv, fnd1, "-nxy"); /* view xvgr file */
- do_view(oenv, fnd2, "-nxy"); /* view xvgr file */
-
- return 0;
-}
*/
#include "legacymodules.h"
+#include <cstdio>
+
#include "gromacs/commandline/cmdlinemodule.h"
#include "gromacs/commandline/cmdlinemodulemanager.h"
-#include "gromacs/onlinehelp/helpwritercontext.h"
-#include "gromacs/utility/exceptions.h"
#include "gromacs/gmxana/gmx_ana.h"
namespace
{
+/*! \brief
+ * Command line module that provides information about obsolescence.
+ *
+ * Prints a message directing the user to a wiki page describing replacement
+ * options.
+ */
+class ObsoleteToolModule : public gmx::CommandLineModuleInterface
+{
+ public:
+ //! Creates an obsolete tool module for a tool with the given name.
+ explicit ObsoleteToolModule(const char *name)
+ : name_(name)
+ {
+ }
+
+ virtual const char *name() const
+ {
+ return name_;
+ }
+ virtual const char *shortDescription() const
+ {
+ return NULL;
+ }
+
+ virtual int run(int /*argc*/, char * /*argv*/[])
+ {
+ printMessage();
+ return 0;
+ }
+ virtual void writeHelp(const gmx::CommandLineHelpContext & /*context*/) const
+ {
+ printMessage();
+ }
+
+ private:
+ void printMessage() const
+ {
+ std::fprintf(stderr,
+ "This tool has been removed from Gromacs 5.0. Please see\n"
+ " http://www.gromacs.org/Documentation/How-tos/Tool_Changes_for_5.0\n"
+ "for ideas how to perform the same tasks with the "
+ "new tools.\n");
+ }
+
+ const char *name_;
+
+};
+
/*! \brief
* Convenience function for creating and registering a module.
*
manager->addModuleCMain(name, shortDescription, mainFunction);
}
+/*! \brief
+ * Convenience function for registering a module for an obsolete tool.
+ *
+ * \param[in] manager Module manager to which to register the module.
+ * \param[in] name Name for the obsolete tool.
+ */
+void registerObsoleteTool(gmx::CommandLineModuleManager *manager,
+ const char *name)
+{
+ gmx::CommandLineModulePointer module(new ObsoleteToolModule(name));
+ manager->addModule(move(module));
+}
+
} // namespace
void registerLegacyModules(gmx::CommandLineModuleManager *manager)
"Calculate distributions and correlations for angles and dihedrals");
registerModule(manager, &gmx_bar, "bar",
"Calculate free energy difference estimates through Bennett's acceptance ratio");
- registerModule(manager, &gmx_bond, "bond",
- "Calculate distances between atoms and bond length distributions");
+ registerObsoleteTool(manager, "bond");
+ registerObsoleteTool(manager, "dist");
+ registerObsoleteTool(manager, "sgangle");
+
registerModule(manager, &gmx_bundle, "bundle",
"Analyze bundles of axes, e.g., helices");
registerModule(manager, &gmx_chi, "chi",
"Compute the total dipole plus fluctuations");
registerModule(manager, &gmx_disre, "disre",
"Analyze distance restraints");
- registerModule(manager, &gmx_dist, "dist",
- "Calculate distances between centers of mass of two groups");
registerModule(manager, &gmx_dos, "dos",
"Analyze density of states and properties based on that");
registerModule(manager, &gmx_dyecoupl, "dyecoupl",
"Compute solvent accessible surface area");
registerModule(manager, &gmx_saxs, "saxs",
"Calculates SAXS structure factors based on Cromer's method");
- registerModule(manager, &gmx_sgangle, "sgangle",
- "Compute the angle and distance between two groups");
registerModule(manager, &gmx_sham, "sham",
"Compute free energies or other histograms from histograms");
registerModule(manager, &gmx_sigeps, "sigeps",
"View a trajectory on an X-Windows terminal");
// TODO: Also include binaries from other directories than src/tools/:
- // "g_xrama|Show animated Ramachandran plots");
// "mdrun|finds a potential energy minimum and calculates the Hessian");
// "mdrun|with -rerun (re)calculates energies for trajectory frames");
}
biod.pnpi.spb.ru / alexxy / gromacs.git / commitdiff