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