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-2008, The GROMACS development team,
6 * check out http://www.gromacs.org for more information.
7 * Copyright (c) 2012,2013, by the GROMACS development team, led by
8 * David van der Spoel, Berk Hess, Erik Lindahl, and including many
9 * others, as listed in the AUTHORS file in the top-level source
10 * directory and at http://www.gromacs.org.
12 * GROMACS is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU Lesser General Public License
14 * as published by the Free Software Foundation; either version 2.1
15 * of the License, or (at your option) any later version.
17 * GROMACS is distributed in the hope that it will be useful,
18 * but WITHOUT ANY WARRANTY; without even the implied warranty of
19 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
20 * Lesser General Public License for more details.
22 * You should have received a copy of the GNU Lesser General Public
23 * License along with GROMACS; if not, see
24 * http://www.gnu.org/licenses, or write to the Free Software Foundation,
25 * Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301 USA.
27 * If you want to redistribute modifications to GROMACS, please
28 * consider that scientific software is very special. Version
29 * control is crucial - bugs must be traceable. We will be happy to
30 * consider code for inclusion in the official distribution, but
31 * derived work must not be called official GROMACS. Details are found
32 * in the README & COPYING files - if they are missing, get the
33 * official version at http://www.gromacs.org.
35 * To help us fund GROMACS development, we humbly ask that you cite
36 * the research papers on the package. Check out http://www.gromacs.org.
50 #include "gmx_fatal.h"
61 t_histo *init_histo(int np,real minx,real maxx)
70 gmx_fatal(FARGS,"minx (%f) should be less than maxx (%f) in init_histo",minx,maxx);
73 h->dx_1 = np/(maxx-minx);
74 h->dx = (maxx-minx)/np;
79 void done_histo(t_histo *h)
86 void add_histo(t_histo *h,real x,real y)
90 n = (x-h->minx)*h->dx_1;
91 if ((n < 0) || (n > h->np))
92 gmx_fatal(FARGS,"Invalid x (%f) in add_histo. Should be in %f - %f",x,h->minx,h->maxx);
97 void dump_histo(t_histo *h,char *fn,char *title,char *xaxis,char *yaxis,
98 int enorm,real norm_fac)
103 for(nn=h->np; (nn > 0); nn--)
106 for(i=0; (i<nn); i++)
109 fp = xvgropen(fn,title,xaxis,yaxis);
110 for( ; (i<nn); i++) {
113 fprintf(fp,"%12f %12f %d\n",h->minx+h->dx*i,h->y[i],h->nh[i]);
116 fprintf(fp,"%12f %12f %d\n",h->minx+h->dx*i,h->y[i]*norm_fac,h->nh[i]);
120 fprintf(fp,"%12f %12f %d\n",
121 h->minx+h->dx*i,h->y[i]*norm_fac/h->nh[i],h->nh[i]);
124 gmx_fatal(FARGS,"Wrong value for enorm (%d)",enorm);
130 /*******************************************************************
132 * Functions to analyse and monitor scattering
134 *******************************************************************/
136 void add_scatter_event(t_ana_scat *scatter,rvec pos,gmx_bool bInel,
139 int np = scatter->np;
141 if (np == scatter->maxp) {
143 srenew(scatter->time,scatter->maxp);
144 srenew(scatter->ekin,scatter->maxp);
145 srenew(scatter->bInel,scatter->maxp);
146 srenew(scatter->pos,scatter->maxp);
148 scatter->time[np] = t;
149 scatter->bInel[np] = np;
150 scatter->ekin[np] = ekin;
151 copy_rvec(pos,scatter->pos[np]);
155 void reset_ana_scat(t_ana_scat *scatter)
160 void done_scatter(t_ana_scat *scatter)
163 sfree(scatter->time);
164 sfree(scatter->ekin);
165 sfree(scatter->bInel);
169 void analyse_scatter(t_ana_scat *scatter,t_histo *hmfp)
174 if (scatter->np > 1) {
175 for(i=1; (i<scatter->np); i++) {
176 rvec_sub(scatter->pos[i],scatter->pos[i-1],dx);
177 add_histo(hmfp,scatter->ekin[i],norm(dx));
182 /*******************************************************************
184 * Functions to analyse structural changes
186 *******************************************************************/
188 t_ana_struct *init_ana_struct(int nstep,int nsave,real timestep,
194 anal->nanal = 1.2*((nstep / nsave)+1);
196 anal->dt = nsave*timestep;
197 snew(anal->t,anal->nanal);
198 snew(anal->maxdist,anal->nanal);
199 snew(anal->d2_com,anal->nanal);
200 snew(anal->d2_origin,anal->nanal);
201 snew(anal->nion,anal->nanal);
204 anal->maxparticle = maxparticle;
207 snew(anal->x[0],maxparticle);
212 void done_ana_struct(t_ana_struct *anal)
217 sfree(anal->maxdist);
219 sfree(anal->d2_origin);
222 for(i=0; (i<anal->nstruct); i++)
227 void reset_ana_struct(t_ana_struct *anal)
231 for(i=0; (i<anal->nanal); i++) {
233 anal->maxdist[i] = 0;
234 clear_rvec(anal->d2_com[i]);
235 clear_rvec(anal->d2_origin[i]);
241 void add_ana_struct(t_ana_struct *total,t_ana_struct *add)
245 if (total->index == 0)
246 total->index = add->index;
247 else if (total->index != add->index)
248 gmx_fatal(FARGS,"Analysis incompatible (total: %d, add: %d) %s, %d",
249 total->index,add->index,__FILE__,__LINE__);
250 for(i=0; (i<total->index); i++) {
251 if (total->t[i] == 0)
252 total->t[i] = add->t[i];
253 else if (total->t[i] != add->t[i])
254 gmx_fatal(FARGS,"Inconsistent times in analysis (%f-%f) %s, %d",
255 total->t[i],add->t[i],__FILE__,__LINE__);
256 if (add->maxdist[i] > total->maxdist[i])
257 total->maxdist[i] = add->maxdist[i];
258 for(m=0; (m<DIM); m++) {
259 total->d2_com[i][m] += add->d2_com[i][m]/add->nion[i];
260 total->d2_origin[i][m] += add->d2_origin[i][m]/add->nion[i];
262 total->nion[i] += add->nion[i];
266 static void do_add_struct(t_ana_struct *anal,int nparticle,rvec x[])
270 if (nparticle > anal->nparticle) {
271 for(i=0; (i<anal->nstruct); i++) {
272 for(j=anal->nparticle; (j<nparticle); j++)
273 copy_rvec(x[j],anal->x[i][j]);
276 anal->nparticle=nparticle;
277 srenew(anal->x,anal->nstruct+1);
278 snew(anal->x[anal->nstruct],anal->maxparticle);
279 for(j=0; (j<nparticle); j++)
280 copy_rvec(x[j],anal->x[anal->nstruct][j]);
284 void analyse_structure(t_ana_struct *anal,real t,rvec center,
285 rvec x[],int nparticle,real charge[])
292 if (j >= anal->nanal)
293 gmx_fatal(FARGS,"Too many points in analyse_structure");
295 anal->maxdist[j] = 0;
298 for(i=0; (i<nparticle); i++) {
307 for(i=0; (i<nparticle); i++) {
309 rvec_sub(x[i],com,dx);
310 for(m=0; (m<DIM); m++) {
311 anal->d2_com[j][m] += sqr(dx[m]);
312 anal->d2_origin[j][m] += sqr(x[i][m]);
314 dx2 = iprod(x[i],x[i]);
316 if (dx1 > anal->maxdist[j])
317 anal->maxdist[j] = dx1;
321 do_add_struct(anal,nparticle,x);
326 void dump_ana_struct(char *rmax,char *nion,char *gyr_com,char *gyr_origin,
327 t_ana_struct *anal,int nsim)
329 FILE *fp,*gp,*hp,*kp;
332 char *legend[] = { "Rg", "RgX", "RgY", "RgZ" };
334 fp = xvgropen(rmax,"rmax","Time (fs)","r (nm)");
335 gp = xvgropen(nion,"N ion","Time (fs)","N ions");
336 hp = xvgropen(gyr_com,"Radius of gyration wrt. C.O.M.",
337 "Time (fs)","Rg (nm)");
338 xvgr_legend(hp,asize(legend),legend);
339 kp = xvgropen(gyr_origin,"Radius of gyration wrt. Origin",
340 "Time (fs)","Rg (nm)");
341 xvgr_legend(kp,asize(legend),legend);
342 for(i=0; (i<anal->index); i++) {
344 fprintf(fp,"%12g %10.3f\n",t,anal->maxdist[i]);
345 fprintf(gp,"%12g %10.3f\n",t,(1.0*anal->nion[i])/nsim-1);
346 d2 = anal->d2_com[i][XX] + anal->d2_com[i][YY] + anal->d2_com[i][ZZ];
347 fprintf(hp,"%12g %10.3f %10.3f %10.3f %10.3f\n",
349 sqrt(anal->d2_com[i][XX]/nsim),
350 sqrt(anal->d2_com[i][YY]/nsim),
351 sqrt(anal->d2_com[i][ZZ]/nsim));
352 d2 = anal->d2_origin[i][XX] + anal->d2_origin[i][YY] + anal->d2_origin[i][ZZ];
353 fprintf(kp,"%12g %10.3f %10.3f %10.3f %10.3f\n",
355 sqrt(anal->d2_origin[i][XX]/nsim),
356 sqrt(anal->d2_origin[i][YY]/nsim),
357 sqrt(anal->d2_origin[i][ZZ]/nsim));
365 void dump_as_pdb(char *pdb,t_ana_struct *anal)
371 kp = ffopen(pdb,"w");
372 for(i=0; (i<anal->nstruct); i++) {
374 fprintf(kp,"MODEL %d time %g fs\n",i+1,t);
375 for(j=0; (j<anal->nparticle); j++) {
376 fprintf(kp,pdbformat,"ATOM",i+1,(j < anal->nion[i]) ? "O" : "N",
378 anal->x[i][j][XX]/100,
379 anal->x[i][j][YY]/100,
380 anal->x[i][j][ZZ]/100);
383 fprintf(kp,"ENDMDL\n");
389 "Coulomb", "Repulsion", "Potential",
390 "EkHole", "EkElectron", "EkLattice", "Kinetic",
394 void add_ana_ener(t_ana_ener *ae,int nn,real e[])
398 /* First time around we are constantly increasing the array size */
400 if (ae->nx == ae->maxx) {
402 srenew(ae->e,ae->maxx);
404 for(i=0; (i<eNR); i++)
405 ae->e[ae->nx][i] = e[i];
409 for(i=0; (i<eNR); i++)
410 ae->e[nn][i] += e[i];
414 void dump_ana_ener(t_ana_ener *ae,int nsim,real dt,char *edump,
421 fac = 1.0/(nsim*ELECTRONVOLT);
422 fp=xvgropen(edump,"Energies","Time (fs)","E (eV)");
423 xvgr_legend(fp,eNR,enms);
424 fprintf(fp,"@ s%d legend \"Ek/Nelec\"\n",eNR);
425 fprintf(fp,"@ type nxy\n");
426 for(i=0; (i<ae->nx); i++) {
427 fprintf(fp,"%10f",1000.0*dt*i);
428 for(j=0; (j<eNR); j++)
429 fprintf(fp," %8.3f",ae->e[i][j]*fac);
430 fprintf(fp," %8.3f\n",ae->e[i][eELECTRON]/(ELECTRONVOLT*total->nion[i]));