?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Redefine the default boolean type to gmx_bool.
[gromacs.git] / src / contrib / ehanal.c
blob754af7e47746ceac8afd2e48d8c907d7b746cd62
1 /*
2 * $Id: ehanal.c,v 1.9.2.3 2008/02/29 07:02:43 spoel Exp $
3 *
4 * This source code is part of
5 *
6 * G R O M A C S
7 *
8 * GROningen MAchine for Chemical Simulations
9 *
10 * VERSION 3.3.3
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2008, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
15
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
20 *
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
27 *
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
30 *
31 * For more info, check our website at http://www.gromacs.org
32 *
33 * And Hey:
34 * Groningen Machine for Chemical Simulation
35 */
36 #ifdef HAVE_CONFIG_H
37 #include <config.h>
38 #endif
39
40 #include <stdlib.h>
41 #include <math.h>
42 #include <string.h>
43 #include "typedefs.h"
44 #include "smalloc.h"
45 #include "macros.h"
46 #include "copyrite.h"
47 #include "statutil.h"
48 #include "gmx_fatal.h"
49 #include "random.h"
50 #include "pdbio.h"
51 #include "futil.h"
52 #include "physics.h"
53 #include "xvgr.h"
54 #include "vec.h"
55 #include "names.h"
56 #include "ehdata.h"
57 #include "pdbio.h"
58
59 t_histo *init_histo(int np,real minx,real maxx)
60 {
61 t_histo *h;
62
63 snew(h,1);
64 snew(h->y,np+1);
65 snew(h->nh,np+1);
66 h->np = np;
67 if (maxx <= minx)
68 gmx_fatal(FARGS,"minx (%f) should be less than maxx (%f) in init_histo",minx,maxx);
69 h->minx = minx;
70 h->maxx = maxx;
71 h->dx_1 = np/(maxx-minx);
72 h->dx = (maxx-minx)/np;
73
74 return h;
75 }
76
77 void done_histo(t_histo *h)
78 {
79 sfree(h->y);
80 sfree(h->nh);
81 h->np = 0;
82 }
83
84 void add_histo(t_histo *h,real x,real y)
85 {
86 int n;
87
88 n = (x-h->minx)*h->dx_1;
89 if ((n < 0) || (n > h->np))
90 gmx_fatal(FARGS,"Invalid x (%f) in add_histo. Should be in %f - %f",x,h->minx,h->maxx);
91 h->y[n] += y;
92 h->nh[n]++;
93 }
94
95 void dump_histo(t_histo *h,char *fn,char *title,char *xaxis,char *yaxis,
96 int enorm,real norm_fac)
97 {
98 FILE *fp;
99 int i,nn;
100
101 for(nn=h->np; (nn > 0); nn--)
102 if (h->nh[nn] != 0)
103 break;
104 for(i=0; (i<nn); i++)
105 if (h->nh[i] > 0)
106 break;
107 fp = xvgropen(fn,title,xaxis,yaxis);
108 for( ; (i<nn); i++) {
109 switch (enorm) {
110 case enormNO:
111 fprintf(fp,"%12f %12f %d\n",h->minx+h->dx*i,h->y[i],h->nh[i]);
112 break;
113 case enormFAC:
114 fprintf(fp,"%12f %12f %d\n",h->minx+h->dx*i,h->y[i]*norm_fac,h->nh[i]);
115 break;
116 case enormNP:
117 if (h->nh[i] > 0)
118 fprintf(fp,"%12f %12f %d\n",
119 h->minx+h->dx*i,h->y[i]*norm_fac/h->nh[i],h->nh[i]);
120 break;
121 default:
122 gmx_fatal(FARGS,"Wrong value for enorm (%d)",enorm);
123 }
124 }
125 ffclose(fp);
126 }
127
128 /*******************************************************************
129 *
130 * Functions to analyse and monitor scattering
131 *
132 *******************************************************************/
133
134 void add_scatter_event(t_ana_scat *scatter,rvec pos,gmx_bool bInel,
135 real t,real ekin)
136 {
137 int np = scatter->np;
138
139 if (np == scatter->maxp) {
140 scatter->maxp += 32;
141 srenew(scatter->time,scatter->maxp);
142 srenew(scatter->ekin,scatter->maxp);
143 srenew(scatter->bInel,scatter->maxp);
144 srenew(scatter->pos,scatter->maxp);
145 }
146 scatter->time[np] = t;
147 scatter->bInel[np] = np;
148 scatter->ekin[np] = ekin;
149 copy_rvec(pos,scatter->pos[np]);
150 scatter->np++;
151 }
152
153 void reset_ana_scat(t_ana_scat *scatter)
154 {
155 scatter->np = 0;
156 }
157
158 void done_scatter(t_ana_scat *scatter)
159 {
160 scatter->np = 0;
161 sfree(scatter->time);
162 sfree(scatter->ekin);
163 sfree(scatter->bInel);
164 sfree(scatter->pos);
165 }
166
167 void analyse_scatter(t_ana_scat *scatter,t_histo *hmfp)
168 {
169 int i,n;
170 rvec dx;
171
172 if (scatter->np > 1) {
173 for(i=1; (i<scatter->np); i++) {
174 rvec_sub(scatter->pos[i],scatter->pos[i-1],dx);
175 add_histo(hmfp,scatter->ekin[i],norm(dx));
176 }
177 }
178 }
179
180 /*******************************************************************
181 *
182 * Functions to analyse structural changes
183 *
184 *******************************************************************/
185
186 t_ana_struct *init_ana_struct(int nstep,int nsave,real timestep,
187 int maxparticle)
188 {
189 t_ana_struct *anal;
190
191 snew(anal,1);
192 anal->nanal = 1.2*((nstep / nsave)+1);
193 anal->index = 0;
194 anal->dt = nsave*timestep;
195 snew(anal->t,anal->nanal);
196 snew(anal->maxdist,anal->nanal);
197 snew(anal->d2_com,anal->nanal);
198 snew(anal->d2_origin,anal->nanal);
199 snew(anal->nion,anal->nanal);
200 anal->nstruct = 1;
201 anal->nparticle = 1;
202 anal->maxparticle = maxparticle;
203 snew(anal->q,1);
204 snew(anal->x,1);
205 snew(anal->x[0],maxparticle);
206
207 return anal;
208 }
209
210 void done_ana_struct(t_ana_struct *anal)
211 {
212 int i;
213
214 sfree(anal->t);
215 sfree(anal->maxdist);
216 sfree(anal->d2_com);
217 sfree(anal->d2_origin);
218 sfree(anal->nion);
219 sfree(anal->q);
220 for(i=0; (i<anal->nstruct); i++)
221 sfree(anal->x[i]);
222 sfree(anal->x);
223 }
224
225 void reset_ana_struct(t_ana_struct *anal)
226 {
227 int i;
228
229 for(i=0; (i<anal->nanal); i++) {
230 anal->t[i] = 0;
231 anal->maxdist[i] = 0;
232 clear_rvec(anal->d2_com[i]);
233 clear_rvec(anal->d2_origin[i]);
234 anal->nion[i] = 0;
235 }
236 anal->index = 0;
237 }
238
239 void add_ana_struct(t_ana_struct *total,t_ana_struct *add)
240 {
241 int i,m;
242
243 if (total->index == 0)
244 total->index = add->index;
245 else if (total->index != add->index)
246 gmx_fatal(FARGS,"Analysis incompatible (total: %d, add: %d) %s, %d",
247 total->index,add->index,__FILE__,__LINE__);
248 for(i=0; (i<total->index); i++) {
249 if (total->t[i] == 0)
250 total->t[i] = add->t[i];
251 else if (total->t[i] != add->t[i])
252 gmx_fatal(FARGS,"Inconsistent times in analysis (%f-%f) %s, %d",
253 total->t[i],add->t[i],__FILE__,__LINE__);
254 if (add->maxdist[i] > total->maxdist[i])
255 total->maxdist[i] = add->maxdist[i];
256 for(m=0; (m<DIM); m++) {
257 total->d2_com[i][m] += add->d2_com[i][m]/add->nion[i];
258 total->d2_origin[i][m] += add->d2_origin[i][m]/add->nion[i];
259 }
260 total->nion[i] += add->nion[i];
261 }
262 }
263
264 static void do_add_struct(t_ana_struct *anal,int nparticle,rvec x[])
265 {
266 int i,j;
267
268 if (nparticle > anal->nparticle) {
269 for(i=0; (i<anal->nstruct); i++) {
270 for(j=anal->nparticle; (j<nparticle); j++)
271 copy_rvec(x[j],anal->x[i][j]);
272 }
273 }
274 anal->nparticle=nparticle;
275 srenew(anal->x,anal->nstruct+1);
276 snew(anal->x[anal->nstruct],anal->maxparticle);
277 for(j=0; (j<nparticle); j++)
278 copy_rvec(x[j],anal->x[anal->nstruct][j]);
279 anal->nstruct++;
280 }
281
282 void analyse_structure(t_ana_struct *anal,real t,rvec center,
283 rvec x[],int nparticle,real charge[])
284 {
285 int i,j,m,nel,n=0;
286 rvec dx,com;
287 real dx2,dx1;
288
289 j = anal->index;
290 if (j >= anal->nanal)
291 gmx_fatal(FARGS,"Too many points in analyse_structure");
292 anal->t[j] = t;
293 anal->maxdist[j] = 0;
294 clear_rvec(com);
295 nel = 0;
296 for(i=0; (i<nparticle); i++) {
297 if (charge[i] < 0) {
298 rvec_inc(com,x[i]);
299 nel++;
300 }
301 }
302 if (nel > 0)
303 for(m=0; (m<3); m++)
304 com[m] /= nel;
305 for(i=0; (i<nparticle); i++) {
306 if (charge[i] < 0) {
307 rvec_sub(x[i],com,dx);
308 for(m=0; (m<DIM); m++) {
309 anal->d2_com[j][m] += sqr(dx[m]);
310 anal->d2_origin[j][m] += sqr(x[i][m]);
311 }
312 dx2 = iprod(x[i],x[i]);
313 dx1 = sqrt(dx2);
314 if (dx1 > anal->maxdist[j])
315 anal->maxdist[j] = dx1;
316 n++;
317 }
318 }
319 do_add_struct(anal,nparticle,x);
320 anal->nion[j] = n;
321 anal->index++;
322 }
323
324 void dump_ana_struct(char *rmax,char *nion,char *gyr_com,char *gyr_origin,
325 t_ana_struct *anal,int nsim)
326 {
327 FILE *fp,*gp,*hp,*kp;
328 int i,j;
329 real t,d2;
330 char *legend[] = { "Rg", "RgX", "RgY", "RgZ" };
331
332 fp = xvgropen(rmax,"rmax","Time (fs)","r (nm)");
333 gp = xvgropen(nion,"N ion","Time (fs)","N ions");
334 hp = xvgropen(gyr_com,"Radius of gyration wrt. C.O.M.",
335 "Time (fs)","Rg (nm)");
336 xvgr_legend(hp,asize(legend),legend);
337 kp = xvgropen(gyr_origin,"Radius of gyration wrt. Origin",
338 "Time (fs)","Rg (nm)");
339 xvgr_legend(kp,asize(legend),legend);
340 for(i=0; (i<anal->index); i++) {
341 t = 1000*anal->t[i];
342 fprintf(fp,"%12g %10.3f\n",t,anal->maxdist[i]);
343 fprintf(gp,"%12g %10.3f\n",t,(1.0*anal->nion[i])/nsim-1);
344 d2 = anal->d2_com[i][XX] + anal->d2_com[i][YY] + anal->d2_com[i][ZZ];
345 fprintf(hp,"%12g %10.3f %10.3f %10.3f %10.3f\n",
346 t,sqrt(d2/nsim),
347 sqrt(anal->d2_com[i][XX]/nsim),
348 sqrt(anal->d2_com[i][YY]/nsim),
349 sqrt(anal->d2_com[i][ZZ]/nsim));
350 d2 = anal->d2_origin[i][XX] + anal->d2_origin[i][YY] + anal->d2_origin[i][ZZ];
351 fprintf(kp,"%12g %10.3f %10.3f %10.3f %10.3f\n",
352 t,sqrt(d2/nsim),
353 sqrt(anal->d2_origin[i][XX]/nsim),
354 sqrt(anal->d2_origin[i][YY]/nsim),
355 sqrt(anal->d2_origin[i][ZZ]/nsim));
356 }
357 ffclose(hp);
358 ffclose(gp);
359 ffclose(fp);
360 ffclose(kp);
361 }
362
363 void dump_as_pdb(char *pdb,t_ana_struct *anal)
364 {
365 FILE *kp;
366 int i,j;
367 real t;
368
369 kp = ffopen(pdb,"w");
370 for(i=0; (i<anal->nstruct); i++) {
371 t = 1000*anal->t[i];
372 fprintf(kp,"MODEL %d time %g fs\n",i+1,t);
373 for(j=0; (j<anal->nparticle); j++) {
374 fprintf(kp,pdbformat,"ATOM",i+1,(j < anal->nion[i]) ? "O" : "N",
375 "PLS",' ',1,
376 anal->x[i][j][XX]/100,
377 anal->x[i][j][YY]/100,
378 anal->x[i][j][ZZ]/100);
379 fprintf(kp,"\n");
380 }
381 fprintf(kp,"ENDMDL\n");
382 }
383 ffclose(kp);
384 }
385
386 char *enms[eNR] = {
387 "Coulomb", "Repulsion", "Potential",
388 "EkHole", "EkElectron", "EkLattice", "Kinetic",
389 "Total"
390 };
391
392 void add_ana_ener(t_ana_ener *ae,int nn,real e[])
393 {
394 int i;
395
396 /* First time around we are constantly increasing the array size */
397 if (nn >= ae->nx) {
398 if (ae->nx == ae->maxx) {
399 ae->maxx += 1024;
400 srenew(ae->e,ae->maxx);
401 }
402 for(i=0; (i<eNR); i++)
403 ae->e[ae->nx][i] = e[i];
404 ae->nx++;
405 }
406 else {
407 for(i=0; (i<eNR); i++)
408 ae->e[nn][i] += e[i];
409 }
410 }
411
412 void dump_ana_ener(t_ana_ener *ae,int nsim,real dt,char *edump,
413 t_ana_struct *total)
414 {
415 FILE *fp;
416 int i,j;
417 real fac;
418
419 fac = 1.0/(nsim*ELECTRONVOLT);
420 fp=xvgropen(edump,"Energies","Time (fs)","E (eV)");
421 xvgr_legend(fp,eNR,enms);
422 fprintf(fp,"@ s%d legend \"Ek/Nelec\"\n",eNR);
423 fprintf(fp,"@ type nxy\n");
424 for(i=0; (i<ae->nx); i++) {
425 fprintf(fp,"%10f",1000.0*dt*i);
426 for(j=0; (j<eNR); j++)
427 fprintf(fp," %8.3f",ae->e[i][j]*fac);
428 fprintf(fp," %8.3f\n",ae->e[i][eELECTRON]/(ELECTRONVOLT*total->nion[i]));
429 }
430 fprintf(fp,"&\n");
431 ffclose(fp);
432 }
433
The ultimate molecular dynamics simulation package