?  search:
summary | log | graphiclog1 | graphiclog2 | commit | commitdiff | tree | refs | edit | fork
blame | history | raw | HEAD
Rename ffopen and ffclose to gmx_ff*.
[gromacs/AngularHB.git] / src / contrib / pmetest.c
blob63d6ef14ff3d6a23ea2d6c78d4e30845508805ed
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.99_development_20071104
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-2006, 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 "typedefs.h"
40 #include "macros.h"
41 #include "smalloc.h"
42 #include "copyrite.h"
43 #include "main.h"
44 #include "nrnb.h"
45 #include "txtdump.h"
46 #include "gromacs/fileio/tpxio.h"
47 #include "gromacs/commandline/pargs.h"
48 #include "gromacs/fileio/futil.h"
49 #include "gmx_fatal.h"
50 #include "vec.h"
51 #include "mdatoms.h"
52 #include "coulomb.h"
53 #include "nsb.h"
54 #include "rmpbc.h"
55 #include "pme.h"
56 #include "force.h"
57 #include "xvgr.h"
58 #include "pbc.h"
59
60 #include "gromacs/utility/gmxmpi.h"
61
62 #include "block_tx.h"
63
64 rvec *xptr=NULL;
65
66 static int comp_xptr(const void *a,const void *b)
67 {
68 int va,vb;
69 real dx;
70
71 va = *(int *)a;
72 vb = *(int *)b;
73
74 if ((dx = (xptr[va][XX] - xptr[vb][XX])) < 0)
75 return -1;
76 else if (dx > 0)
77 return 1;
78 else
79 return 0;
80 }
81
82 static void do_my_pme(FILE *fp,real tm,gmx_bool bVerbose,t_inputrec *ir,
83 rvec x[],rvec xbuf[],rvec f[],
84 real charge[],real qbuf[],real qqbuf[],
85 matrix box,gmx_bool bSort,
86 t_commrec *cr,t_nsborder *nsb,t_nrnb *nrnb,
87 t_block *excl,real qtot,
88 t_forcerec *fr,int index[],FILE *fp_xvg,
89 int ngroups,unsigned short cENER[])
90 {
91 real ener,vcorr,q,xx,dvdl=0,vdip,vcharge;
92 tensor vir,vir_corr,vir_tot;
93 rvec mu_tot[2];
94 int i,m,ii,ig,jg;
95 real **epme,*qptr;
96
97 /* Initiate local variables */
98 fr->f_el_recip = f;
99 clear_mat(vir);
100 clear_mat(vir_corr);
101
102 if (ngroups > 1) {
103 fprintf(fp,"There are %d energy groups\n",ngroups);
104 snew(epme,ngroups);
105 for(i=0; (i<ngroups); i++)
106 snew(epme[i],ngroups);
107 }
108
109 /* Put x is in the box, this part needs to be parallellized properly */
110 /*put_atoms_in_box(box,nsb->natoms,x);*/
111 /* Here sorting of X (and q) is done.
112 * Alternatively, one could just put the atoms in one of the
113 * cr->nnodes slabs. That is much cheaper than sorting.
114 */
115 for(i=0; (i<nsb->natoms); i++)
116 index[i] = i;
117 if (bSort) {
118 xptr = x;
119 qsort(index,nsb->natoms,sizeof(index[0]),comp_xptr);
120 xptr = NULL; /* To trap unintentional use of the ptr */
121 }
122 /* After sorting we only need the part that is to be computed on
123 * this processor. We also compute the mu_tot here (system dipole)
124 */
125 clear_rvec(mu_tot[0]);
126 for(i=START(nsb); (i<START(nsb)+HOMENR(nsb)); i++) {
127 ii = index[i];
128 q = charge[ii];
129 qbuf[i] = q;
130 for(m=0; (m<DIM); m++) {
131 xx = x[ii][m];
132 xbuf[i][m] = xx;
133 mu_tot[0][m] += q*xx;
134 }
135 clear_rvec(f[ii]);
136 }
137 copy_rvec(mu_tot[0],mu_tot[1]);
138 if (debug) {
139 pr_rvec(debug,0,"qbuf",qbuf,nsb->natoms,TRUE);
140 pr_rvecs(debug,0,"xbuf",xbuf,nsb->natoms);
141 pr_rvecs(debug,0,"box",box,DIM);
142 }
143 for(ig=0; (ig<ngroups); ig++) {
144 for(jg=ig; (jg<ngroups); jg++) {
145 if (ngroups > 1) {
146 for(i=START(nsb); (i<START(nsb)+HOMENR(nsb)); i++) {
147 if ((cENER[i] == ig) || (cENER[i] == jg))
148 qqbuf[i] = qbuf[i];
149 else
150 qqbuf[i] = 0;
151 }
152 qptr = qqbuf;
153 }
154 else
155 qptr = qbuf;
156 ener = do_pme(fp,bVerbose,ir,xbuf,f,qptr,qptr,box,cr,
157 nsb,nrnb,vir,fr->ewaldcoeff,FALSE,0,&dvdl,FALSE);
158 vcorr = ewald_LRcorrection(fp,nsb,cr,fr,qptr,qptr,excl,xbuf,box,mu_tot,
159 ir->ewald_geometry,ir->epsilon_surface,
160 0,&dvdl,&vdip,&vcharge);
161 gmx_sum(1,&ener,cr);
162 gmx_sum(1,&vcorr,cr);
163 if (ngroups > 1)
164 epme[ig][jg] = ener+vcorr;
165 }
166 }
167 if (ngroups > 1) {
168 if (fp_xvg)
169 fprintf(fp_xvg,"%10.3f",tm);
170 for(ig=0; (ig<ngroups); ig++) {
171 for(jg=ig; (jg<ngroups); jg++) {
172 if (ig != jg)
173 epme[ig][jg] -= epme[ig][ig]+epme[jg][jg];
174 if (fp_xvg)
175 fprintf(fp_xvg," %12.5e",epme[ig][jg]);
176 }
177 }
178 if (fp_xvg)
179 fprintf(fp_xvg,"\n");
180 }
181 else {
182 fprintf(fp,"Time: %10.3f Energy: %12.5e Correction: %12.5e Total: %12.5e\n",
183 tm,ener,vcorr,ener+vcorr);
184 if (fp_xvg)
185 fprintf(fp_xvg,"%10.3f %12.5e %12.5e %12.5e\n",tm,ener+vcorr,vdip,vcharge);
186 if (bVerbose) {
187 m_add(vir,vir_corr,vir_tot);
188 gmx_sum(9,vir_tot[0],cr);
189 pr_rvecs(fp,0,"virial",vir_tot,DIM);
190 }
191 fflush(fp);
192 }
193 }
194
195 int main(int argc,char *argv[])
196 {
197 static char *desc[] = {
198 "The [TT]pmetest[tt] program tests the scaling of the PME code. When only given",
199 "a [TT].tpr[tt] file it will compute PME for one frame. When given a trajectory",
200 "it will do so for all the frames in the trajectory. Before the PME",
201 "routine is called the coordinates are sorted along the X-axis.[PAR]",
202 "As an extra service to the public the program can also compute",
203 "long-range Coulomb energies for components of the system. When the",
204 "[TT]-groups[tt] flag is given to the program the energy groups",
205 "from the [TT].tpr[tt] file will be read, and half an energy matrix computed."
206 };
207 t_commrec *cr,*mcr;
208 static t_filenm fnm[] = {
209 { efTPX, NULL, NULL, ffREAD },
210 { efTRN, "-o", NULL, ffWRITE },
211 { efLOG, "-g", "pme", ffWRITE },
212 { efTRX, "-f", NULL, ffOPTRD },
213 { efXVG, "-x", "ener-pme", ffWRITE }
214 };
215 #define NFILE asize(fnm)
216
217 /* Command line options ! */
218 static gmx_bool bVerbose=FALSE;
219 static gmx_bool bOptFFT=FALSE;
220 static gmx_bool bSort=FALSE;
221 static int ewald_geometry=eewg3D;
222 static int nnodes=1;
223 static int pme_order=0;
224 static rvec grid = { -1, -1, -1 };
225 static real rc = 0.0;
226 static real dtol = 0.0;
227 static gmx_bool bGroups = FALSE;
228 static t_pargs pa[] = {
229 { "-np", FALSE, etINT, {&nnodes},
230 "Number of nodes, must be the same as used for [TT]grompp[tt]" },
231 { "-v", FALSE, etBOOL,{&bVerbose},
232 "Be loud and noisy" },
233 { "-sort", FALSE, etBOOL,{&bSort},
234 "Sort coordinates. Crucial for domain decomposition." },
235 { "-grid", FALSE, etRVEC,{&grid},
236 "Number of grid cells in X, Y, Z dimension (if -1 use from [TT].tpr[tt])" },
237 { "-order", FALSE, etINT, {&pme_order},
238 "Order of the PME spreading algorithm" },
239 { "-groups", FALSE, etBOOL, {&bGroups},
240 "Compute half an energy matrix based on the energy groups in your [TT].tpr[tt] file" },
241 { "-rc", FALSE, etREAL, {&rc},
242 "Rcoulomb for Ewald summation" },
243 { "-tol", FALSE, etREAL, {&dtol},
244 "Tolerance for Ewald summation" }
245 };
246 FILE *fp;
247 t_inputrec *ir;
248 t_topology top;
249 t_tpxheader tpx;
250 t_nrnb nrnb;
251 t_nsborder *nsb;
252 t_forcerec *fr;
253 t_mdatoms *mdatoms;
254 char title[STRLEN];
255 int natoms,step,status,i,ncg,root;
256 real t,lambda,ewaldcoeff,qtot;
257 rvec *x,*f,*xbuf;
258 int *index;
259 gmx_bool bCont;
260 real *charge,*qbuf,*qqbuf;
261 matrix box;
262
263 /* Start the actual parallel code if necessary */
264 cr = init_par(&argc,&argv);
265 root = 0;
266
267 if (MASTER(cr))
268 CopyRight(stderr,argv[0]);
269
270 /* Parse command line on all processors, arguments are passed on in
271 * init_par (see above)
272 */
273 parse_common_args(&argc,argv,
274 PCA_NOEXIT_ON_ARGS | PCA_BE_NICE |
275 PCA_CAN_SET_DEFFNM | (MASTER(cr) ? 0 : PCA_QUIET),
276 NFILE,fnm,asize(pa),pa,asize(desc),desc,0,NULL);
277
278 #ifndef GMX_MPI
279 if (nnodes > 1)
280 gmx_fatal(FARGS,"GROMACS compiled without MPI support - can't do parallel runs");
281 #endif
282
283 /* Open log files on all processors */
284 open_log(ftp2fn(efLOG,NFILE,fnm),cr);
285 snew(ir,1);
286
287 if (MASTER(cr)) {
288 /* Read tpr file etc. */
289 read_tpxheader(ftp2fn(efTPX,NFILE,fnm),&tpx,FALSE,NULL,NULL);
290 snew(x,tpx.natoms);
291 read_tpx(ftp2fn(efTPX,NFILE,fnm),&step,&t,&lambda,ir,
292 box,&natoms,x,NULL,NULL,&top);
293 /* Charges */
294 qtot = 0;
295 snew(charge,natoms);
296 for(i=0; (i<natoms); i++) {
297 charge[i] = top.atoms.atom[i].q;
298 qtot += charge[i];
299 }
300
301 /* Grid stuff */
302 if (opt2parg_bSet("-grid",asize(pa),pa)) {
303 ir->nkx = grid[XX];
304 ir->nky = grid[YY];
305 ir->nkz = grid[ZZ];
306 }
307 /* Check command line parameters for consistency */
308 if ((ir->nkx <= 0) || (ir->nky <= 0) || (ir->nkz <= 0))
309 gmx_fatal(FARGS,"PME grid = %d %d %d",ir->nkx,ir->nky,ir->nkz);
310 if (opt2parg_bSet("-rc",asize(pa),pa))
311 ir->rcoulomb = rc;
312 if (ir->rcoulomb <= 0)
313 gmx_fatal(FARGS,"rcoulomb should be > 0 (not %f)",ir->rcoulomb);
314 if (opt2parg_bSet("-order",asize(pa),pa))
315 ir->pme_order = pme_order;
316 if (ir->pme_order <= 0)
317 gmx_fatal(FARGS,"pme_order should be > 0 (not %d)",ir->pme_order);
318 if (opt2parg_bSet("-tol",asize(pa),pa))
319 ir->ewald_rtol = dtol;
320 if (ir->ewald_rtol <= 0)
321 gmx_fatal(FARGS,"ewald_tol should be > 0 (not %f)",ir->ewald_rtol);
322 }
323 else {
324 init_top(&top);
325 }
326
327 /* Add parallellization code here */
328 snew(nsb,1);
329 if (MASTER(cr)) {
330 ncg = top.blocks[ebCGS].multinr[0];
331 for(i=0; (i<cr->nnodes-1); i++)
332 top.blocks[ebCGS].multinr[i] = min(ncg,(ncg*(i+1))/cr->nnodes);
333 for( ; (i<MAXNODES); i++)
334 top.blocks[ebCGS].multinr[i] = ncg;
335 }
336 if (PAR(cr)) {
337 /* Set some variables to zero to avoid core dumps */
338 ir->opts.ngtc = ir->opts.ngacc = ir->opts.ngfrz = ir->opts.ngener = 0;
339 #ifdef GMX_MPI
340 /* Distribute the data over processors */
341 MPI_Bcast(&natoms,1,MPI_INT,root,MPI_COMM_WORLD);
342 MPI_Bcast(ir,sizeof(*ir),MPI_BYTE,root,MPI_COMM_WORLD);
343 MPI_Bcast(&qtot,1,GMX_MPI_REAL,root,MPI_COMM_WORLD);
344 #endif
345
346 /* Call some dedicated communication routines, master sends n-1 times */
347 if (MASTER(cr)) {
348 for(i=1; (i<cr->nnodes); i++) {
349 mv_block(i,&(top.blocks[ebCGS]));
350 mv_block(i,&(top.atoms.excl));
351 }
352 }
353 else {
354 ld_block(root,&(top.blocks[ebCGS]));
355 ld_block(root,&(top.atoms.excl));
356 }
357 if (!MASTER(cr)) {
358 snew(charge,natoms);
359 snew(x,natoms);
360 }
361 #ifdef GMX_MPI
362 MPI_Bcast(charge,natoms,GMX_MPI_REAL,root,MPI_COMM_WORLD);
363 #endif
364 }
365 ewaldcoeff = calc_ewaldcoeff(ir->rcoulomb,ir->ewald_rtol);
366
367
368 if (bVerbose)
369 pr_inputrec(stdlog,0,"Inputrec",ir);
370
371 /* Allocate memory for temp arrays etc. */
372 snew(xbuf,natoms);
373 snew(f,natoms);
374 snew(qbuf,natoms);
375 snew(qqbuf,natoms);
376 snew(index,natoms);
377
378 /* Initialize the PME code */
379 init_pme(stdlog,cr,ir->nkx,ir->nky,ir->nkz,ir->pme_order,
380 natoms,FALSE,bOptFFT,ewald_geometry);
381
382 /* MFlops accounting */
383 init_nrnb(&nrnb);
384
385 /* Initialize the work division */
386 calc_nsb(stdlog,&(top.blocks[ebCGS]),cr->nnodes,nsb,0);
387 nsb->nodeid = cr->nodeid;
388 print_nsb(stdlog,"pmetest",nsb);
389
390 /* Initiate forcerec */
391 mdatoms = atoms2md(stdlog,&top.atoms,ir->opts.nFreeze,ir->eI,
392 ir->delta_t,0,ir->opts.tau_t,FALSE,FALSE);
393 snew(fr,1);
394 init_forcerec(stdlog,fr,ir,&top,cr,mdatoms,nsb,box,FALSE,NULL,NULL,FALSE);
395
396 /* First do PME based on coordinates in tpr file, send them to
397 * other processors if needed.
398 */
399 if (MASTER(cr))
400 fprintf(stdlog,"-----\n"
401 "Results based on tpr file %s\n",ftp2fn(efTPX,NFILE,fnm));
402 #ifdef GMX_MPI
403 if (PAR(cr)) {
404 MPI_Bcast(x[0],natoms*DIM,GMX_MPI_REAL,root,MPI_COMM_WORLD);
405 MPI_Bcast(box[0],DIM*DIM,GMX_MPI_REAL,root,MPI_COMM_WORLD);
406 MPI_Bcast(&t,1,GMX_MPI_REAL,root,MPI_COMM_WORLD);
407 }
408 #endif
409 do_my_pme(stdlog,0,bVerbose,ir,x,xbuf,f,charge,qbuf,qqbuf,box,bSort,
410 cr,nsb,&nrnb,&(top.atoms.excl),qtot,fr,index,NULL,
411 bGroups ? ir->opts.ngener : 1,mdatoms->cENER);
412
413 /* If we have a trajectry file, we will read the frames in it and compute
414 * the PME energy.
415 */
416 if (ftp2bSet(efTRX,NFILE,fnm)) {
417 fprintf(stdlog,"-----\n"
418 "Results based on trx file %s\n",ftp2fn(efTRX,NFILE,fnm));
419 if (MASTER(cr)) {
420 sfree(x);
421 natoms = read_first_x(&status,ftp2fn(efTRX,NFILE,fnm),&t,&x,box);
422 if (natoms != top.atoms.nr)
423 gmx_fatal(FARGS,"natoms in trx = %d, in tpr = %d",natoms,top.atoms.nr);
424 fp = xvgropen(ftp2fn(efXVG,NFILE,fnm),"PME Energy","Time (ps)","E (kJ/mol)");
425 }
426 else
427 fp = NULL;
428 do {
429 /* Send coordinates, box and time to the other nodes */
430 #ifdef GMX_MPI
431 if (PAR(cr)) {
432 MPI_Bcast(x[0],natoms*DIM,GMX_MPI_REAL,root,MPI_COMM_WORLD);
433 MPI_Bcast(box[0],DIM*DIM,GMX_MPI_REAL,root,MPI_COMM_WORLD);
434 MPI_Bcast(&t,1,GMX_MPI_REAL,root,MPI_COMM_WORLD);
435 }
436 #endif
437 rm_pbc(&top.idef,nsb->natoms,box,x,x);
438 /* Call the PME wrapper function */
439 do_my_pme(stdlog,t,bVerbose,ir,x,xbuf,f,charge,qbuf,qqbuf,box,bSort,cr,
440 nsb,&nrnb,&(top.atoms.excl),qtot,fr,index,fp,
441 bGroups ? ir->opts.ngener : 1,mdatoms->cENER);
442 /* Only the master processor reads more data */
443 if (MASTER(cr))
444 bCont = read_next_x(status,&t,natoms,x,box);
445 /* Check whether we need to continue */
446 #ifdef GMX_MPI
447 if (PAR(cr))
448 MPI_Bcast(&bCont,1,MPI_INT,root,MPI_COMM_WORLD);
449 #endif
450
451 } while (bCont);
452
453 /* Finish I/O, close files */
454 if (MASTER(cr)) {
455 close_trx(status);
456 gmx_ffclose(fp);
457 }
458 }
459
460 if (bVerbose) {
461 /* Do some final I/O about performance, might be useful in debugging */
462 fprintf(stdlog,"-----\n");
463 print_nrnb(stdlog,&nrnb);
464 }
465
466 /* Finish the parallel stuff */
467 if (gmx_parallel_env_initialized())
468 gmx_finalize(cr);
469
470 /* Thank the audience, as usual */
471 if (MASTER(cr))
472 gmx_thanx(stderr);
473
474 return 0;
475 }
476