?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Merge branch 'master' of git://git.gromacs.org/gromacs
[gromacs/adressmacs.git] / src / tools / gmx_pme_error.c
blobc0396989d6c7eda030075d16e0acaaa66fe15ff8
1 /* $Id: gmx_tune_pme.c 9 2009-08-11 09:43:30Z dommert $
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 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
10 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
11 * Copyright (c) 2001-2008, The GROMACS development team,
12 * check out http://www.gromacs.org for more information.
13
14 * This program is free software; you can redistribute it and/or
15 * modify it under the terms of the GNU General Public License
16 * as published by the Free Software Foundation; either version 2
17 * of the License, or (at your option) any later version.
18 *
19 * If you want to redistribute modifications, please consider that
20 * scientific software is very special. Version control is crucial -
21 * bugs must be traceable. We will be happy to consider code for
22 * inclusion in the official distribution, but derived work must not
23 * be called official GROMACS. Details are found in the README & COPYING
24 * files - if they are missing, get the official version at www.gromacs.org.
25 *
26 * To help us fund GROMACS development, we humbly ask that you cite
27 * the papers on the package - you can find them in the top README file.
28 *
29 * For more info, check our website at http://www.gromacs.org
30 *
31 * And Hey:
32 * Gallium Rubidium Oxygen Manganese Argon Carbon Silicon
33 */
34 #include "statutil.h"
35 #include "typedefs.h"
36 #include "smalloc.h"
37 #include "vec.h"
38 #include "copyrite.h"
39 #include "tpxio.h"
40 #include "string2.h"
41 #include "readinp.h"
42 #include "calcgrid.h"
43 #include "checkpoint.h"
44 #include "gmx_ana.h"
45 #include "gmx_random.h"
46 #include "physics.h"
47 #include "mdatoms.h"
48 #include "coulomb.h"
49 #include "mtop_util.h"
50 #include "network.h"
51 #include "main.h"
52
53 /* We use the same defines as in mvdata.c here */
54 #define block_bc(cr, d) gmx_bcast( sizeof(d), &(d),(cr))
55 #define nblock_bc(cr,nr,d) gmx_bcast((nr)*sizeof((d)[0]), (d),(cr))
56 #define snew_bc(cr,d,nr) { if (!MASTER(cr)) snew((d),(nr)); }
57 /* #define TAKETIME */
58 /* #define DEBUG */
59 enum {
60 ddnoSEL, ddnoINTERLEAVE, ddnoPP_PME, ddnoCARTESIAN, ddnoNR
61 };
62
63 /* Enum for situations that can occur during log file parsing */
64 enum {
65 eParselogOK,
66 eParselogNotFound,
67 eParselogNoPerfData,
68 eParselogTerm,
69 eParselogResetProblem,
70 eParselogNr
71 };
72
73
74 typedef struct
75 {
76 int nPMEnodes; /* number of PME only nodes used in this test */
77 int nx, ny, nz; /* DD grid */
78 int guessPME; /* if nPMEnodes == -1, this is the guessed number of PME nodes */
79 float *Gcycles; /* This can contain more than one value if doing multiple tests */
80 float Gcycles_Av;
81 float *ns_per_day;
82 float ns_per_day_Av;
83 float *PME_f_load; /* PME mesh/force load average*/
84 float PME_f_load_Av; /* Average average ;) ... */
85 char *mdrun_cmd_line; /* Mdrun command line used for this test */
86 } t_perf;
87
88
89 typedef struct
90 {
91 gmx_large_int_t orig_sim_steps; /* Number of steps to be done in the real simulation */
92 int n_entries; /* Number of entries in arrays */
93 real volume; /* The volume of the box */
94 matrix recipbox; /* The reciprocal box */
95 int natoms; /* The number of atoms in the MD system */
96 real *fac; /* The scaling factor */
97 real *rcoulomb; /* The coulomb radii [0...nr_inputfiles] */
98 real *rvdw; /* The vdW radii */
99 int *nkx, *nky, *nkz; /* Number of k vectors in each spatial dimension */
100 real *fourier_sp; /* Fourierspacing */
101 real *ewald_rtol; /* Real space tolerance for Ewald, determines */
102 /* the real/reciprocal space relative weight */
103 real *ewald_beta; /* Splitting parameter [1/nm] */
104 real fracself; /* fraction of particles for SI error */
105 real q2all; /* sum ( q ^2 ) */
106 real q2allnr; /* nr of charges */
107 int *pme_order; /* Interpolation order for PME (bsplines) */
108 char **fn_out; /* Name of the output tpr file */
109 real *e_dir; /* Direct space part of PME error with these settings */
110 real *e_rec; /* Reciprocal space part of PME error */
111 gmx_bool bTUNE; /* flag for tuning */
112 } t_inputinfo;
113
114
115 /* Returns TRUE when atom is charged */
116 static gmx_bool is_charge(real charge)
117 {
118 if (charge*charge > GMX_REAL_EPS)
119 return TRUE;
120 else
121 return FALSE;
122 }
123
124
125 /* calculate charge density */
126 static void calc_q2all(
127 gmx_mtop_t *mtop, /* molecular topology */
128 real *q2all, real *q2allnr)
129 {
130 int imol,iatom; /* indices for loops */
131 real q2_all=0; /* Sum of squared charges */
132 int nrq_mol; /* Number of charges in a single molecule */
133 int nrq_all; /* Total number of charges in the MD system */
134 real nrq_all_r; /* No of charges in real format */
135 real qi,q2_mol;
136 gmx_moltype_t *molecule;
137 gmx_molblock_t *molblock;
138
139 #ifdef DEBUG
140 fprintf(stderr, "\nCharge density:\n");
141 #endif
142 q2_all = 0.0; /* total q squared */
143 nrq_all = 0; /* total number of charges in the system */
144 for (imol=0; imol<mtop->nmolblock; imol++) /* Loop over molecule types */
145 {
146 q2_mol=0.0; /* q squared value of this molecule */
147 nrq_mol=0; /* number of charges this molecule carries */
148 molecule = &(mtop->moltype[imol]);
149 molblock = &(mtop->molblock[imol]);
150 for (iatom=0; iatom<molblock->natoms_mol; iatom++) /* Loop over atoms in this molecule */
151 {
152 qi = molecule->atoms.atom[iatom].q; /* Charge of this atom */
153 /* Is this charge worth to be considered? */
154 if (is_charge(qi))
155 {
156 q2_mol += qi*qi;
157 nrq_mol++;
158 }
159 }
160 /* Multiply with the number of molecules present of this type and add */
161 q2_all += q2_mol*molblock->nmol;
162 nrq_all += nrq_mol*molblock->nmol;
163 #ifdef DEBUG
164 fprintf(stderr, "Molecule %2d (%5d atoms) q2_mol=%10.3e nr.mol.charges=%5d (%6dx) q2_all=%10.3e tot.charges=%d\n",
165 imol,molblock->natoms_mol,q2_mol,nrq_mol,molblock->nmol,q2_all,nrq_all);
166 #endif
167 }
168 nrq_all_r = nrq_all;
169
170 *q2all=q2_all;
171 *q2allnr=nrq_all;
172
173 }
174
175
176 /* Estimate the direct space part error of the SPME Ewald sum */
177 static real estimate_direct(
178 t_inputinfo *info
179 )
180 {
181 real e_dir=0; /* Error estimate */
182 real beta=0; /* Splitting parameter (1/nm) */
183 real r_coulomb=0; /* Cut-off in direct space */
184
185
186 beta = info->ewald_beta[0];
187 r_coulomb = info->rcoulomb[0];
188
189 e_dir = 2.0 * info->q2all * gmx_invsqrt( info->q2allnr * r_coulomb * info->volume );
190 e_dir *= exp (-beta*beta*r_coulomb*r_coulomb);
191
192 return ONE_4PI_EPS0*e_dir;
193 }
194
195 #define SUMORDER 6
196
197 /* the following 4 functions determine polynomials required for the reciprocal error estimate */
198
199 static inline real eps_poly1(
200 real m, /* grid coordinate in certain direction */
201 real K, /* grid size in corresponding direction */
202 real n) /* spline interpolation order of the SPME */
203 {
204 int i;
205 real nom=0; /* nominator */
206 real denom=0; /* denominator */
207 real tmp=0;
208
209 if ( m == 0.0 )
210 return 0.0 ;
211
212 for(i=-SUMORDER ; i<0 ; i++)
213 {
214 tmp=m / K + i;
215 tmp*=2.0*M_PI;
216 nom+=pow( tmp , -n );
217 }
218
219 for(i=SUMORDER ; i>0 ; i--)
220 {
221 tmp=m / K + i;
222 tmp*=2.0*M_PI;
223 nom+=pow( tmp , -n );
224 }
225
226 tmp=m / K;
227 tmp*=2.0*M_PI;
228 denom=pow( tmp , -n )+nom;
229
230 return -nom/denom;
231
232 }
233
234 static inline real eps_poly2(
235 real m, /* grid coordinate in certain direction */
236 real K, /* grid size in corresponding direction */
237 real n) /* spline interpolation order of the SPME */
238 {
239 int i;
240 real nom=0; /* nominator */
241 real denom=0; /* denominator */
242 real tmp=0;
243
244 if ( m == 0.0 )
245 return 0.0 ;
246
247 for(i=-SUMORDER ; i<0 ; i++)
248 {
249 tmp=m / K + i;
250 tmp*=2.0*M_PI;
251 nom+=pow( tmp , -2.0*n );
252 }
253
254 for(i=SUMORDER ; i>0 ; i--)
255 {
256 tmp=m / K + i;
257 tmp*=2.0*M_PI;
258 nom+=pow( tmp , -2.0*n );
259 }
260
261 for(i=-SUMORDER ; i<SUMORDER+1 ; i++)
262 {
263 tmp=m / K + i;
264 tmp*=2.0*M_PI;
265 denom+=pow( tmp , -n );
266 }
267 tmp=eps_poly1(m,K,n);
268 return nom / denom / denom + tmp*tmp ;
269
270 }
271
272 static inline real eps_poly3(
273 real m, /* grid coordinate in certain direction */
274 real K, /* grid size in corresponding direction */
275 real n) /* spline interpolation order of the SPME */
276 {
277 int i;
278 real nom=0; /* nominator */
279 real denom=0; /* denominator */
280 real tmp=0;
281
282 if ( m == 0.0 )
283 return 0.0 ;
284
285 for(i=-SUMORDER ; i<0 ; i++)
286 {
287 tmp=m / K + i;
288 tmp*=2.0*M_PI;
289 nom+= i * pow( tmp , -2.0*n );
290 }
291
292 for(i=SUMORDER ; i>0 ; i--)
293 {
294 tmp=m / K + i;
295 tmp*=2.0*M_PI;
296 nom+= i * pow( tmp , -2.0*n );
297 }
298
299 for(i=-SUMORDER ; i<SUMORDER+1 ; i++)
300 {
301 tmp=m / K + i;
302 tmp*=2.0*M_PI;
303 denom+=pow( tmp , -n );
304 }
305
306 return 2.0 * M_PI * nom / denom / denom;
307
308 }
309
310 static inline real eps_poly4(
311 real m, /* grid coordinate in certain direction */
312 real K, /* grid size in corresponding direction */
313 real n) /* spline interpolation order of the SPME */
314 {
315 int i;
316 real nom=0; /* nominator */
317 real denom=0; /* denominator */
318 real tmp=0;
319
320 if ( m == 0.0 )
321 return 0.0 ;
322
323 for(i=-SUMORDER ; i<0 ; i++)
324 {
325 tmp=m / K + i;
326 tmp*=2.0*M_PI;
327 nom+= i * i * pow( tmp , -2.0*n );
328 }
329
330 for(i=SUMORDER ; i>0 ; i--)
331 {
332 tmp=m / K + i;
333 tmp*=2.0*M_PI;
334 nom+= i * i * pow( tmp , -2.0*n );
335 }
336
337 for(i=-SUMORDER ; i<SUMORDER+1 ; i++)
338 {
339 tmp=m / K + i;
340 tmp*=2.0*M_PI;
341 denom+=pow( tmp , -n );
342 }
343
344 return 4.0 * M_PI * M_PI * nom / denom / denom;
345
346 }
347
348 static inline real eps_self(
349 real m, /* grid coordinate in certain direction */
350 real K, /* grid size in corresponding direction */
351 rvec rboxv, /* reciprocal box vector */
352 real n, /* spline interpolation order of the SPME */
353 rvec x) /* coordinate of charge */
354 {
355 int i;
356 real tmp=0; /* temporary variables for computations */
357 real tmp1=0; /* temporary variables for computations */
358 real tmp2=0; /* temporary variables for computations */
359 real rcoord=0; /* coordinate in certain reciprocal space direction */
360 real nom=0; /* nominator */
361 real denom=0; /* denominator */
362
363
364 if ( m == 0.0 )
365 return 0.0 ;
366
367 rcoord=iprod(rboxv,x);
368
369
370 for(i=-SUMORDER;i<0;i++)
371 {
372 tmp=-sin(2.0 * M_PI * i * K * rcoord);
373 tmp1=2.0 * M_PI * m / K + 2.0 * M_PI * i;
374 tmp2=pow(tmp1,-1.0*n);
375 nom+=tmp * tmp2 * i;
376 denom+=tmp2;
377 }
378
379 for(i=SUMORDER;i>0;i--)
380 {
381 tmp=-sin(2.0 * M_PI * i * K * rcoord);
382 tmp1=2.0 * M_PI * m / K + 2.0 * M_PI * i;
383 tmp2=pow(tmp1,-1.0*n);
384 nom+=tmp * tmp2 * i;
385 denom+=tmp2;
386 }
387
388
389 tmp=2.0 * M_PI * m / K;
390 tmp1=pow(tmp,-1.0*n);
391 denom+=tmp1;
392
393 return 2.0 * M_PI * nom / denom * K ;
394
395 }
396
397 #undef SUMORDER
398
399 /* The following routine is just a copy from pme.c */
400
401 static void calc_recipbox(matrix box,matrix recipbox)
402 {
403 /* Save some time by assuming upper right part is zero */
404
405 real tmp=1.0/(box[XX][XX]*box[YY][YY]*box[ZZ][ZZ]);
406
407 recipbox[XX][XX]=box[YY][YY]*box[ZZ][ZZ]*tmp;
408 recipbox[XX][YY]=0;
409 recipbox[XX][ZZ]=0;
410 recipbox[YY][XX]=-box[YY][XX]*box[ZZ][ZZ]*tmp;
411 recipbox[YY][YY]=box[XX][XX]*box[ZZ][ZZ]*tmp;
412 recipbox[YY][ZZ]=0;
413 recipbox[ZZ][XX]=(box[YY][XX]*box[ZZ][YY]-box[YY][YY]*box[ZZ][XX])*tmp;
414 recipbox[ZZ][YY]=-box[ZZ][YY]*box[XX][XX]*tmp;
415 recipbox[ZZ][ZZ]=box[XX][XX]*box[YY][YY]*tmp;
416 }
417
418
419 /* Estimate the reciprocal space part error of the SPME Ewald sum. */
420 static real estimate_reciprocal(
421 t_inputinfo *info,
422 rvec x[], /* array of particles */
423 real q[], /* array of charges */
424 int nr, /* number of charges = size of the charge array */
425 FILE *fp_out,
426 gmx_bool bVerbose,
427 unsigned int seed, /* The seed for the random number generator */
428 int *nsamples, /* Return the number of samples used if Monte Carlo
429 * algorithm is used for self energy error estimate */
430 t_commrec *cr)
431 {
432 real e_rec=0; /* reciprocal error estimate */
433 real e_rec1=0; /* Error estimate term 1*/
434 real e_rec2=0; /* Error estimate term 2*/
435 real e_rec3=0; /* Error estimate term 3 */
436 real e_rec3x=0; /* part of Error estimate term 3 in x */
437 real e_rec3y=0; /* part of Error estimate term 3 in y */
438 real e_rec3z=0; /* part of Error estimate term 3 in z */
439 int i,ci;
440 int nx,ny,nz; /* grid coordinates */
441 real q2_all=0; /* sum of squared charges */
442 rvec gridpx; /* reciprocal grid point in x direction*/
443 rvec gridpxy; /* reciprocal grid point in x and y direction*/
444 rvec gridp; /* complete reciprocal grid point in 3 directions*/
445 rvec tmpvec; /* template to create points from basis vectors */
446 rvec tmpvec2; /* template to create points from basis vectors */
447 real coeff=0; /* variable to compute coefficients of the error estimate */
448 real coeff2=0; /* variable to compute coefficients of the error estimate */
449 real tmp=0; /* variables to compute different factors from vectors */
450 real tmp1=0;
451 real tmp2=0;
452 gmx_bool bFraction;
453
454 /* Random number generator */
455 gmx_rng_t rng=NULL;
456 int *numbers=NULL;
457
458 /* Index variables for parallel work distribution */
459 int startglobal,stopglobal;
460 int startlocal, stoplocal;
461 int x_per_core;
462 int xtot;
463
464 #ifdef TAKETIME
465 double t0=0.0;
466 double t1=0.0;
467 #endif
468
469 rng=gmx_rng_init(seed);
470
471 clear_rvec(gridpx);
472 clear_rvec(gridpxy);
473 clear_rvec(gridp);
474 clear_rvec(tmpvec);
475 clear_rvec(tmpvec2);
476
477 for(i=0;i<nr;i++)
478 {
479 q2_all += q[i]*q[i];
480 }
481
482 /* Calculate indices for work distribution */
483 startglobal=-info->nkx[0]/2;
484 stopglobal = info->nkx[0]/2;
485 xtot = stopglobal*2+1;
486 if (PAR(cr))
487 {
488 x_per_core = ceil((real)xtot / (real)cr->nnodes);
489 startlocal = startglobal + x_per_core*cr->nodeid;
490 stoplocal = startlocal + x_per_core -1;
491 if (stoplocal > stopglobal)
492 stoplocal = stopglobal;
493 }
494 else
495 {
496 startlocal = startglobal;
497 stoplocal = stopglobal;
498 x_per_core = xtot;
499 }
500 /*
501 #ifdef GMX_MPI
502 MPI_Barrier(MPI_COMM_WORLD);
503 #endif
504 */
505
506 #ifdef GMX_MPI
507 #ifdef TAKETIME
508 if (MASTER(cr))
509 t0 = MPI_Wtime();
510 #endif
511 #endif
512
513 if (MASTER(cr)){
514
515 fprintf(stderr, "Calculating reciprocal error part 1 ...");
516
517 }
518
519 for(nx=startlocal; nx<=stoplocal; nx++)
520 {
521 svmul(nx,info->recipbox[XX],gridpx);
522 for(ny=-info->nky[0]/2; ny<info->nky[0]/2+1; ny++)
523 {
524 svmul(ny,info->recipbox[YY],tmpvec);
525 rvec_add(gridpx,tmpvec,gridpxy);
526 for(nz=-info->nkz[0]/2; nz<info->nkz[0]/2+1; nz++)
527 {
528 if ( 0 == nx && 0 == ny && 0 == nz )
529 continue;
530 svmul(nz,info->recipbox[ZZ],tmpvec);
531 rvec_add(gridpxy,tmpvec,gridp);
532 tmp=norm2(gridp);
533 coeff=exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] ) ;
534 coeff/= 2.0 * M_PI * info->volume * tmp;
535 coeff2=tmp ;
536
537
538 tmp=eps_poly2(nx,info->nkx[0],info->pme_order[0]);
539 tmp+=eps_poly2(ny,info->nkx[0],info->pme_order[0]);
540 tmp+=eps_poly2(nz,info->nkx[0],info->pme_order[0]);
541
542 tmp1=eps_poly1(nx,info->nkx[0],info->pme_order[0]);
543 tmp2=eps_poly1(ny,info->nky[0],info->pme_order[0]);
544
545 tmp+=2.0 * tmp1 * tmp2;
546
547 tmp1=eps_poly1(nz,info->nkz[0],info->pme_order[0]);
548 tmp2=eps_poly1(ny,info->nky[0],info->pme_order[0]);
549
550 tmp+=2.0 * tmp1 * tmp2;
551
552 tmp1=eps_poly1(nz,info->nkz[0],info->pme_order[0]);
553 tmp2=eps_poly1(nx,info->nkx[0],info->pme_order[0]);
554
555 tmp+=2.0 * tmp1 * tmp2;
556
557 tmp1=eps_poly1(nx,info->nkx[0],info->pme_order[0]);
558 tmp1+=eps_poly1(ny,info->nky[0],info->pme_order[0]);
559 tmp1+=eps_poly1(nz,info->nkz[0],info->pme_order[0]);
560
561 tmp+= tmp1 * tmp1;
562
563 e_rec1+= 32.0 * M_PI * M_PI * coeff * coeff * coeff2 * tmp * q2_all * q2_all / nr ;
564
565 tmp1=eps_poly3(nx,info->nkx[0],info->pme_order[0]);
566 tmp1*=info->nkx[0];
567 tmp2=iprod(gridp,info->recipbox[XX]);
568
569 tmp=tmp1*tmp2;
570
571 tmp1=eps_poly3(ny,info->nky[0],info->pme_order[0]);
572 tmp1*=info->nky[0];
573 tmp2=iprod(gridp,info->recipbox[YY]);
574
575 tmp+=tmp1*tmp2;
576
577 tmp1=eps_poly3(nz,info->nkz[0],info->pme_order[0]);
578 tmp1*=info->nkz[0];
579 tmp2=iprod(gridp,info->recipbox[ZZ]);
580
581 tmp+=tmp1*tmp2;
582
583 tmp*=4.0 * M_PI;
584
585 tmp1=eps_poly4(nx,info->nkx[0],info->pme_order[0]);
586 tmp1*=norm2(info->recipbox[XX]);
587 tmp1*=info->nkx[0] * info->nkx[0];
588
589 tmp+=tmp1;
590
591 tmp1=eps_poly4(ny,info->nky[0],info->pme_order[0]);
592 tmp1*=norm2(info->recipbox[YY]);
593 tmp1*=info->nky[0] * info->nky[0];
594
595 tmp+=tmp1;
596
597 tmp1=eps_poly4(nz,info->nkz[0],info->pme_order[0]);
598 tmp1*=norm2(info->recipbox[ZZ]);
599 tmp1*=info->nkz[0] * info->nkz[0];
600
601 tmp+=tmp1;
602
603 e_rec2+= 4.0 * coeff * coeff * tmp * q2_all * q2_all / nr ;
604
605 }
606 }
607 if (MASTER(cr))
608 fprintf(stderr, "\rCalculating reciprocal error part 1 ... %3.0f%%", 100.0*(nx-startlocal+1)/(x_per_core));
609
610 }
611
612 if (MASTER(cr))
613 fprintf(stderr, "\n");
614
615 /* Use just a fraction of all charges to estimate the self energy error term? */
616 bFraction = (info->fracself > 0.0) && (info->fracself < 1.0);
617
618 if (bFraction)
619 {
620 /* Here xtot is the number of samples taken for the Monte Carlo calculation
621 * of the average of term IV of equation 35 in Wang2010. Round up to a
622 * number of samples that is divisible by the number of nodes */
623 x_per_core = ceil(info->fracself * nr / (real)cr->nnodes);
624 xtot = x_per_core * cr->nnodes;
625 }
626 else
627 {
628 /* In this case we use all nr particle positions */
629 xtot = nr;
630 x_per_core = ceil( (real)xtot / (real)cr->nnodes );
631 }
632
633 startlocal = x_per_core * cr->nodeid;
634 stoplocal = min(startlocal + x_per_core, xtot); /* min needed if xtot == nr */
635
636 if (bFraction)
637 {
638 /* Make shure we get identical results in serial and parallel. Therefore,
639 * take the sample indices from a single, global random number array that
640 * is constructed on the master node and that only depends on the seed */
641 snew(numbers, xtot);
642 if (MASTER(cr))
643 {
644 for (i=0; i<xtot; i++)
645 {
646 numbers[i] = floor(gmx_rng_uniform_real(rng) * nr );
647 }
648 }
649 /* Broadcast the random number array to the other nodes */
650 if (PAR(cr))
651 {
652 nblock_bc(cr,xtot,numbers);
653 }
654
655 if (bVerbose && MASTER(cr))
656 {
657 fprintf(stdout, "Using %d sample%s to approximate the self interaction error term",
658 xtot, xtot==1?"":"s");
659 if (PAR(cr))
660 fprintf(stdout, " (%d sample%s per node)", x_per_core, x_per_core==1?"":"s");
661 fprintf(stdout, ".\n");
662 }
663 }
664
665 /* Return the number of positions used for the Monte Carlo algorithm */
666 *nsamples = xtot;
667
668 for(i=startlocal;i<stoplocal;i++)
669 {
670 e_rec3x=0;
671 e_rec3y=0;
672 e_rec3z=0;
673
674 if (bFraction)
675 {
676 /* Randomly pick a charge */
677 ci = numbers[i];
678 }
679 else
680 {
681 /* Use all charges */
682 ci = i;
683 }
684
685 /* for(nx=startlocal; nx<=stoplocal; nx++)*/
686 for(nx=-info->nkx[0]/2; nx<info->nkx[0]/2+1; nx++)
687 {
688 svmul(nx,info->recipbox[XX],gridpx);
689 for(ny=-info->nky[0]/2; ny<info->nky[0]/2+1; ny++)
690 {
691 svmul(ny,info->recipbox[YY],tmpvec);
692 rvec_add(gridpx,tmpvec,gridpxy);
693 for(nz=-info->nkz[0]/2; nz<info->nkz[0]/2+1; nz++)
694 {
695
696 if ( 0 == nx && 0 == ny && 0 == nz)
697 continue;
698
699 svmul(nz,info->recipbox[ZZ],tmpvec);
700 rvec_add(gridpxy,tmpvec,gridp);
701 tmp=norm2(gridp);
702 coeff=exp(-1.0 * M_PI * M_PI * tmp / info->ewald_beta[0] / info->ewald_beta[0] );
703 coeff/= tmp ;
704 e_rec3x+=coeff*eps_self(nx,info->nkx[0],info->recipbox[XX],info->pme_order[0],x[ci]);
705 e_rec3y+=coeff*eps_self(ny,info->nky[0],info->recipbox[YY],info->pme_order[0],x[ci]);
706 e_rec3z+=coeff*eps_self(nz,info->nkz[0],info->recipbox[ZZ],info->pme_order[0],x[ci]);
707
708 }
709 }
710 }
711
712 clear_rvec(tmpvec2);
713
714 svmul(e_rec3x,info->recipbox[XX],tmpvec);
715 rvec_inc(tmpvec2,tmpvec);
716 svmul(e_rec3y,info->recipbox[YY],tmpvec);
717 rvec_inc(tmpvec2,tmpvec);
718 svmul(e_rec3z,info->recipbox[ZZ],tmpvec);
719 rvec_inc(tmpvec2,tmpvec);
720
721 e_rec3 += q[ci]*q[ci]*q[ci]*q[ci]*norm2(tmpvec2) / ( xtot * M_PI * info->volume * M_PI * info->volume);
722 if (MASTER(cr)){
723 fprintf(stderr, "\rCalculating reciprocal error part 2 ... %3.0f%%",
724 100.0*(i+1)/stoplocal);
725
726 }
727 }
728
729 if (MASTER(cr))
730 fprintf(stderr, "\n");
731
732
733 #ifdef TAKETIME
734 if (MASTER(cr))
735 {
736 t1= MPI_Wtime() - t0;
737 fprintf(fp_out, "Recip. err. est. took : %lf s\n", t1);
738 }
739 #endif
740
741 #ifdef DEBUG
742 if (PAR(cr))
743 {
744 fprintf(stderr, "Node %3d: nx=[%3d...%3d] e_rec3=%e\n",
745 cr->nodeid, startlocal, stoplocal, e_rec3);
746 }
747 #endif
748
749 if (PAR(cr))
750 {
751 gmx_sum(1,&e_rec1,cr);
752 gmx_sum(1,&e_rec2,cr);
753 gmx_sum(1,&e_rec3,cr);
754 }
755
756 /* e_rec1*=8.0 * q2_all / info->volume / info->volume / nr ;
757 e_rec2*= q2_all / M_PI / M_PI / info->volume / info->volume / nr ;
758 e_rec3/= M_PI * M_PI * info->volume * info->volume * nr ;
759 */
760 e_rec=sqrt(e_rec1+e_rec2+e_rec3);
761
762
763 return ONE_4PI_EPS0 * e_rec;
764 }
765
766
767 /* Allocate memory for the inputinfo struct: */
768 static void create_info(t_inputinfo *info)
769 {
770 snew(info->fac , info->n_entries);
771 snew(info->rcoulomb , info->n_entries);
772 snew(info->rvdw , info->n_entries);
773 snew(info->nkx , info->n_entries);
774 snew(info->nky , info->n_entries);
775 snew(info->nkz , info->n_entries);
776 snew(info->fourier_sp, info->n_entries);
777 snew(info->ewald_rtol, info->n_entries);
778 snew(info->ewald_beta, info->n_entries);
779 snew(info->pme_order , info->n_entries);
780 snew(info->fn_out , info->n_entries);
781 snew(info->e_dir , info->n_entries);
782 snew(info->e_rec , info->n_entries);
783 }
784
785
786 /* Allocate and fill an array with coordinates and charges,
787 * returns the number of charges found
788 */
789 static int prepare_x_q(real *q[], rvec *x[], gmx_mtop_t *mtop, rvec x_orig[], t_commrec *cr)
790 {
791 int i,anr_global;
792 int nq; /* number of charged particles */
793 t_atom *atom;
794
795
796 if (MASTER(cr))
797 {
798 snew(*q, mtop->natoms);
799 snew(*x, mtop->natoms);
800 nq=0;
801 for (i=0; i<mtop->natoms; i++)
802 {
803 anr_global = i;
804 gmx_mtop_atomnr_to_atom(mtop,anr_global,&atom);
805 if (is_charge(atom->q))
806 {
807 (*q)[nq] = atom->q;
808 (*x)[nq][XX] = x_orig[i][XX];
809 (*x)[nq][YY] = x_orig[i][YY];
810 (*x)[nq][ZZ] = x_orig[i][ZZ];
811 nq++;
812 }
813 }
814 /* Give back some unneeded memory */
815 srenew(*q, nq);
816 srenew(*x, nq);
817 }
818 /* Broadcast x and q in the parallel case */
819 if (PAR(cr))
820 {
821 /* Transfer the number of charges */
822 block_bc(cr,nq);
823 snew_bc(cr, *x, nq);
824 snew_bc(cr, *q, nq);
825 nblock_bc(cr,nq,*x);
826 nblock_bc(cr,nq,*q);
827 }
828
829 return nq;
830 }
831
832
833
834 /* Read in the tpr file and save information we need later in info */
835 static void read_tpr_file(const char *fn_sim_tpr, t_inputinfo *info, t_state *state, gmx_mtop_t *mtop, t_inputrec *ir, real user_beta, real fracself)
836 {
837 read_tpx_state(fn_sim_tpr,ir,state,NULL,mtop);
838
839 /* The values of the original tpr input file are save in the first
840 * place [0] of the arrays */
841 info->orig_sim_steps = ir->nsteps;
842 info->pme_order[0] = ir->pme_order;
843 info->rcoulomb[0] = ir->rcoulomb;
844 info->rvdw[0] = ir->rvdw;
845 info->nkx[0] = ir->nkx;
846 info->nky[0] = ir->nky;
847 info->nkz[0] = ir->nkz;
848 info->ewald_rtol[0] = ir->ewald_rtol;
849 info->fracself = fracself;
850 if (user_beta > 0)
851 info->ewald_beta[0] = user_beta;
852 else
853 info->ewald_beta[0] = calc_ewaldcoeff(info->rcoulomb[0],info->ewald_rtol[0]);
854
855 /* Check if PME was chosen */
856 if (EEL_PME(ir->coulombtype) == FALSE)
857 gmx_fatal(FARGS, "Can only do optimizations for simulations with PME");
858
859 /* Check if rcoulomb == rlist, which is necessary for PME */
860 if (!(ir->rcoulomb == ir->rlist))
861 gmx_fatal(FARGS, "PME requires rcoulomb (%f) to be equal to rlist (%f).", ir->rcoulomb, ir->rlist);
862 }
863
864
865 /* Transfer what we need for parallelizing the reciprocal error estimate */
866 static void bcast_info(t_inputinfo *info, t_commrec *cr)
867 {
868 nblock_bc(cr, info->n_entries, info->nkx);
869 nblock_bc(cr, info->n_entries, info->nky);
870 nblock_bc(cr, info->n_entries, info->nkz);
871 nblock_bc(cr, info->n_entries, info->ewald_beta);
872 nblock_bc(cr, info->n_entries, info->pme_order);
873 nblock_bc(cr, info->n_entries, info->e_dir);
874 nblock_bc(cr, info->n_entries, info->e_rec);
875 block_bc(cr, info->volume);
876 block_bc(cr, info->recipbox);
877 block_bc(cr, info->natoms);
878 block_bc(cr, info->fracself);
879 block_bc(cr, info->bTUNE);
880 block_bc(cr, info->q2all);
881 block_bc(cr, info->q2allnr);
882 }
883
884
885 /* Estimate the error of the SPME Ewald sum. This estimate is based upon
886 * a) a homogeneous distribution of the charges
887 * b) a total charge of zero.
888 */
889 static void estimate_PME_error(t_inputinfo *info, t_state *state,
890 gmx_mtop_t *mtop, FILE *fp_out, gmx_bool bVerbose, unsigned int seed,
891 t_commrec *cr)
892 {
893 rvec *x=NULL; /* The coordinates */
894 real *q=NULL; /* The charges */
895 real edir=0.0; /* real space error */
896 real erec=0.0; /* reciprocal space error */
897 real derr=0.0; /* difference of real and reciprocal space error */
898 real derr0=0.0; /* difference of real and reciprocal space error */
899 real beta=0.0; /* splitting parameter beta */
900 real beta0=0.0; /* splitting parameter beta */
901 int ncharges; /* The number of atoms with charges */
902 int nsamples; /* The number of samples used for the calculation of the
903 * self-energy error term */
904 int i=0;
905
906 if (MASTER(cr))
907 fprintf(fp_out, "\n--- PME ERROR ESTIMATE ---\n");
908
909 /* Prepare an x and q array with only the charged atoms */
910 ncharges = prepare_x_q(&q, &x, mtop, state->x, cr);
911 if (MASTER(cr))
912 {
913 calc_q2all(mtop, &(info->q2all), &(info->q2allnr));
914 info->ewald_rtol[0]=gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
915 /* Write some info to log file */
916 fprintf(fp_out, "Box volume : %g nm^3\n", info->volume);
917 fprintf(fp_out, "Number of charged atoms : %d (total atoms %d)\n",ncharges, info->natoms);
918 fprintf(fp_out, "Coulomb radius : %g nm\n", info->rcoulomb[0]);
919 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
920 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
921 fprintf(fp_out, "Interpolation order : %d\n", info->pme_order[0]);
922 fprintf(fp_out, "Fourier grid (nx,ny,nz) : %d x %d x %d\n",
923 info->nkx[0],info->nky[0],info->nkz[0]);
924 fflush(fp_out);
925
926 }
927
928 if (PAR(cr))
929 bcast_info(info, cr);
930
931
932 /* Calculate direct space error */
933 info->e_dir[0] = estimate_direct(info);
934
935 /* Calculate reciprocal space error */
936 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
937 seed, &nsamples, cr);
938
939 if (PAR(cr))
940 bcast_info(info, cr);
941
942 if (MASTER(cr))
943 {
944 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
945 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
946 fprintf(fp_out, "Self-energy error term was estimated using %d samples\n", nsamples);
947 fflush(fp_out);
948 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
949 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
950 }
951
952 i=0;
953
954 if (info->bTUNE)
955 {
956 if(MASTER(cr))
957 fprintf(stderr,"Starting tuning ...\n");
958 edir=info->e_dir[0];
959 erec=info->e_rec[0];
960 derr0=edir-erec;
961 beta0=info->ewald_beta[0];
962 if (derr>0.0)
963 info->ewald_beta[0]+=0.1;
964 else
965 info->ewald_beta[0]-=0.1;
966 info->e_dir[0] = estimate_direct(info);
967 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
968 seed, &nsamples, cr);
969
970 if (PAR(cr))
971 bcast_info(info, cr);
972
973
974 edir=info->e_dir[0];
975 erec=info->e_rec[0];
976 derr=edir-erec;
977 while ( fabs(derr/min(erec,edir)) > 1e-4)
978 {
979
980 beta=info->ewald_beta[0];
981 beta-=derr*(info->ewald_beta[0]-beta0)/(derr-derr0);
982 beta0=info->ewald_beta[0];
983 info->ewald_beta[0]=beta;
984 derr0=derr;
985
986 info->e_dir[0] = estimate_direct(info);
987 info->e_rec[0] = estimate_reciprocal(info, x, q, ncharges, fp_out, bVerbose,
988 seed, &nsamples, cr);
989
990 if (PAR(cr))
991 bcast_info(info, cr);
992
993 edir=info->e_dir[0];
994 erec=info->e_rec[0];
995 derr=edir-erec;
996
997 if (MASTER(cr))
998 {
999 i++;
1000 fprintf(stderr,"difference between real and rec. space error (step %d): %g\n",i,fabs(derr));
1001 fprintf(stderr,"old beta: %f\n",beta0);
1002 fprintf(stderr,"new beta: %f\n",beta);
1003 }
1004 }
1005
1006 info->ewald_rtol[0]=gmx_erfc(info->rcoulomb[0]*info->ewald_beta[0]);
1007
1008 if (MASTER(cr))
1009 {
1010 /* Write some info to log file */
1011 fflush(fp_out);
1012 fprintf(fp_out, "========= After tuning ========\n");
1013 fprintf(fp_out, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1014 fprintf(fp_out, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1015 fprintf(stderr, "Direct space error est. : %10.3e kJ/(mol*nm)\n", info->e_dir[0]);
1016 fprintf(stderr, "Reciprocal sp. err. est.: %10.3e kJ/(mol*nm)\n", info->e_rec[0]);
1017 fprintf(fp_out, "Ewald_rtol : %g\n", info->ewald_rtol[0]);
1018 fprintf(fp_out, "Ewald parameter beta : %g\n", info->ewald_beta[0]);
1019 fflush(fp_out);
1020
1021 }
1022
1023 }
1024
1025 }
1026
1027
1028 int gmx_pme_error(int argc,char *argv[])
1029 {
1030 const char *desc[] = {
1031 "g_pme_error estimates the error of the electrostatic forces",
1032 "if using the sPME algorithm. The flag [TT]-tune[tt] will determine",
1033 "the splitting parameter such that the error is equally",
1034 "distributed over the real and reciprocal space part.",
1035 "The part of the error that stems from self interaction of the particles "
1036 "is computationally demanding. However, a good a approximation is to",
1037 "just use a fraction of the particles for this term which can be",
1038 "indicated by the flag [TT]-self[tt].[PAR]",
1039 };
1040
1041 real fs=0.0; /* 0 indicates: not set by the user */
1042 real user_beta=-1.0;
1043 real fracself=1.0;
1044 t_inputinfo info;
1045 t_state state; /* The state from the tpr input file */
1046 gmx_mtop_t mtop; /* The topology from the tpr input file */
1047 t_inputrec *ir=NULL; /* The inputrec from the tpr file */
1048 FILE *fp=NULL;
1049 t_commrec *cr;
1050 unsigned long PCA_Flags;
1051 gmx_bool bTUNE=FALSE;
1052 gmx_bool bVerbose=FALSE;
1053 int seed=0;
1054
1055
1056 static t_filenm fnm[] = {
1057 { efTPX, "-s", NULL, ffREAD },
1058 { efOUT, "-o", "error", ffWRITE },
1059 { efTPX, "-so", "tuned", ffOPTWR }
1060 };
1061
1062 output_env_t oenv=NULL;
1063
1064 t_pargs pa[] = {
1065 { "-beta", FALSE, etREAL, {&user_beta},
1066 "If positive, overwrite ewald_beta from tpr file with this value" },
1067 { "-tune", FALSE, etBOOL, {&bTUNE},
1068 "Tune the splitting parameter such that the error is equally distributed between real and reciprocal space" },
1069 { "-self", FALSE, etREAL, {&fracself},
1070 "If between 0.0 and 1.0, determine self interaction error from just this fraction of the charged particles" },
1071 { "-seed", FALSE, etINT, {&seed},
1072 "Random number seed used for Monte Carlo algorithm when -self is set to a value between 0.0 and 1.0" },
1073 { "-v", FALSE, etBOOL, {&bVerbose},
1074 "Be loud and noisy" }
1075 };
1076
1077
1078 #define NFILE asize(fnm)
1079
1080 cr = init_par(&argc,&argv);
1081
1082 #ifdef GMX_MPI
1083 if (PAR(cr))
1084 MPI_Barrier(MPI_COMM_WORLD);
1085 #endif
1086
1087 if (MASTER(cr))
1088 CopyRight(stderr,argv[0]);
1089
1090 PCA_Flags = PCA_NOEXIT_ON_ARGS;
1091 PCA_Flags |= (MASTER(cr) ? 0 : PCA_QUIET);
1092
1093 parse_common_args(&argc,argv,PCA_Flags,
1094 NFILE,fnm,asize(pa),pa,asize(desc),desc,
1095 0,NULL,&oenv);
1096
1097 if (!bTUNE)
1098 bTUNE = opt2bSet("-so",NFILE,fnm);
1099
1100 info.n_entries = 1;
1101
1102 /* Allocate memory for the inputinfo struct: */
1103 create_info(&info);
1104 info.fourier_sp[0] = fs;
1105
1106 /* Read in the tpr file and open logfile for reading */
1107 if (MASTER(cr))
1108 {
1109 snew(ir,1);
1110 read_tpr_file(opt2fn("-s",NFILE,fnm), &info, &state, &mtop, ir, user_beta,fracself);
1111
1112 fp=fopen(opt2fn("-o",NFILE,fnm),"w");
1113 }
1114
1115 /* Check consistency if the user provided fourierspacing */
1116 if (fs > 0 && MASTER(cr))
1117 {
1118 /* Recalculate the grid dimensions using fourierspacing from user input */
1119 info.nkx[0] = 0;
1120 info.nky[0] = 0;
1121 info.nkz[0] = 0;
1122 calc_grid(stdout,state.box,info.fourier_sp[0],&(info.nkx[0]),&(info.nky[0]),&(info.nkz[0]));
1123 if ( (ir->nkx != info.nkx[0]) || (ir->nky != info.nky[0]) || (ir->nkz != info.nkz[0]) )
1124 gmx_fatal(FARGS, "Wrong fourierspacing %f nm, input file grid = %d x %d x %d, computed grid = %d x %d x %d",
1125 fs,ir->nkx,ir->nky,ir->nkz,info.nkx[0],info.nky[0],info.nkz[0]);
1126 }
1127
1128 /* Estimate (S)PME force error */
1129
1130 /* Determine the volume of the simulation box */
1131 if (MASTER(cr))
1132 {
1133 info.volume = det(state.box);
1134 calc_recipbox(state.box,info.recipbox);
1135 info.natoms = mtop.natoms;
1136 info.bTUNE = bTUNE;
1137 }
1138
1139 if (PAR(cr))
1140 bcast_info(&info, cr);
1141
1142 /* Get an error estimate of the input tpr file and do some tuning if requested */
1143 estimate_PME_error(&info, &state, &mtop, fp, bVerbose, seed, cr);
1144
1145 if (MASTER(cr))
1146 {
1147 /* Write out optimized tpr file if requested */
1148 if ( opt2bSet("-so",NFILE,fnm) || bTUNE )
1149 {
1150 ir->ewald_rtol=info.ewald_rtol[0];
1151 write_tpx_state(opt2fn("-so",NFILE,fnm),ir,&state,&mtop);
1152 }
1153 please_cite(fp,"Wang2010");
1154 fclose(fp);
1155 }
1156
1157 if (gmx_parallel_env_initialized())
1158 {
1159 gmx_finalize();
1160 }
1161
1162 return 0;
1163 }
got merge a long time ago