3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * Green Red Orange Magenta Azure Cyan Skyblue
49 #include "gmx_fatal.h"
62 static t_charge *mk_charge(t_atoms *atoms,t_block *cgs,int *nncg)
66 int i,j,ncg,resnr,anr;
69 /* Find the charged groups */
71 for(i=0; (i<cgs->nr); i++) {
73 for(j=cgs->index[i]; (j<cgs->index[i+1]); j++) {
76 if (fabs(qq) > 1.0e-5) {
81 resnr=atoms->atom[anr].resind;
82 sprintf(buf,"%s%d-%d",
83 *(atoms->resinfo[resnr].name),
84 atoms->resinfo[resnr].nr,
86 cg[ncg].label=strdup(buf);
92 for(i=0; (i<ncg); i++) {
93 printf("CG: %10s Q: %6g Atoms:",
95 for(j=cgs->index[cg[i].cg]; (j<cgs->index[cg[i].cg+1]); j++)
103 static real calc_dist(t_pbc *pbc,rvec x[],t_block *cgs,int icg,int jcg)
107 real d2,mindist2=1000;
109 for(i=cgs->index[icg]; (i<cgs->index[icg+1]); i++) {
110 for(j=cgs->index[jcg]; (j<cgs->index[jcg+1]); j++) {
111 pbc_dx(pbc,x[i],x[j],dx);
117 return sqrt(mindist2);
120 int gmx_saltbr(int argc,char *argv[])
122 const char *desc[] = {
123 "[TT]g_saltbr[tt] plots the distance between all combination of charged groups",
124 "as a function of time. The groups are combined in different ways.",
125 "A minimum distance can be given, (ie. a cut-off), then groups",
126 "that are never closer than that distance will not be plotted.[BR]",
127 "Output will be in a number of fixed filenames, [TT]min-min.xvg[tt], [TT]plus-min.xvg[tt]",
128 "and [TT]plus-plus.xvg[tt], or files for every individual ion-pair if the [TT]-sep[tt]",
129 "option is selected. In this case files are named as [TT]sb-ResnameResnr-Atomnr[tt].",
130 "There may be many such files."
132 static gmx_bool bSep=FALSE;
133 static real truncate=1000.0;
135 { "-t", FALSE, etREAL, {&truncate},
137 { "-sep", FALSE, etBOOL, {&bSep},
138 "Use separate files for each interaction (may be MANY)" }
141 { efTRX, "-f", NULL, ffREAD },
142 { efTPX, NULL, NULL, ffREAD },
144 #define NFILE asize(fnm)
147 static const char *title[3] = {
148 "Distance between positively charged groups",
149 "Distance between negatively charged groups",
150 "Distance between oppositely charged groups"
152 static const char *fn[3] = {
163 int i,j,k,m,nnn,teller,ncg,n1,n2,n3,natoms;
176 CopyRight(stderr,argv[0]);
178 parse_common_args(&argc,argv,PCA_CAN_TIME | PCA_BE_NICE,
179 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL,&oenv);
181 top=read_top(ftp2fn(efTPX,NFILE,fnm),&ePBC);
182 cg=mk_charge(&top->atoms,&(top->cgs),&ncg);
185 for(i=0; (i<ncg); i++) {
187 snew(nWithin[i],ncg);
190 natoms=read_first_x(oenv,&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
195 srenew(time,teller+1);
198 set_pbc(&pbc,ePBC,box);
200 for(i=0; (i<ncg); i++) {
201 for(j=i+1; (j<ncg); j++) {
202 srenew(cgdist[i][j],teller+1);
203 cgdist[i][j][teller]=
204 calc_dist(&pbc,x,&(top->cgs),cg[i].cg,cg[j].cg);
205 if (cgdist[i][j][teller] < truncate)
211 } while (read_next_x(oenv,status,&t,natoms,x,box));
212 fprintf(stderr,"\n");
217 for(i=0; (i<ncg); i++)
218 for(j=i+1; (j<ncg); j++) {
220 sprintf(buf,"sb-%s:%s.xvg",cg[i].label,cg[j].label);
221 fp=xvgropen(buf,buf,"Time (ps)","Distance (nm)",oenv);
222 for(k=0; (k<teller); k++)
223 fprintf(fp,"%10g %10g\n",time[k],cgdist[i][j][k]);
232 out[m]=xvgropen(fn[m],title[m],"Time (ps)","Distance (nm)",oenv);
235 for(i=0; (i<ncg); i++) {
237 for(j=i+1; (j<ncg); j++) {
240 sprintf(buf,"%s:%s",cg[i].label,cg[j].label);
249 xvgr_legend(out[nnn],1,(const char**)&buf,oenv);
251 if (output_env_get_xvg_format(oenv) == exvgXMGR) {
252 fprintf(out[nnn],"@ legend string %d \"%s\"\n",nset[nnn],buf);
253 } else if (output_env_get_xvg_format(oenv) == exvgXMGRACE) {
254 fprintf(out[nnn],"@ s%d legend \"%s\"\n",nset[nnn],buf);
262 for(k=0; (k<teller); k++) {
264 fprintf(out[m],"%10g",time[k]);
266 for(i=0; (i<ncg); i++) {
267 for(j=i+1; (j<ncg); j++) {
270 fprintf(out[nnn-1]," %10g",cgdist[i][j][k]);
274 fprintf(out[m],"\n");
276 for(m=0; (m<3); m++) {