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