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 "g_saltbr 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, (eg. the 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, min-min.xvg, plus-min.xvg",
128 "and plus-plus.xvg, 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
++) {