3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
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-2004, The GROMACS development team,
13 * check out http://www.gromacs.org for more information.
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.
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.
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.
30 * For more info, check our website at http://www.gromacs.org
33 * GROningen Mixture of Alchemy and Childrens' Stories
38 #include "gromacs/math/vec.h"
39 #include "gromacs/math/units.h"
42 #include "gromacs/utility/smalloc.h"
43 #include "gromacs/fileio/tpxio.h"
44 #include "gromacs/commandline/pargs.h"
45 #include "gromacs/fileio/writeps.h"
52 #include "gromacs/pbcutil/mshift.h"
56 static real
phi_sr(FILE *log
,int nj
,rvec x
[],real charge
[],real rc
,real r1
,
57 rvec box
, real phi
[],t_block
*excl
,rvec f_sr
[],gmx_bool bOld
)
60 real pp
,r2
,R
,R_1
,R_2
,rc2
;
61 real qi
,qj
,vsr
,eps
,fscal
;
68 for(i
=0; (i
<nj
-1); i
++) {
70 for(j
=i
+1; (j
<nj
); j
++) {
73 for(k
=i1
; (k
<i2
); k
++)
77 r2
=calc_dx2dx(x
[i
],x
[j
],box
,dx
);
84 fscal
= old_f(R
,rc
,r1
)*R_2
;
85 pp
= old_phi(R
,rc
,r1
);
88 fscal
= new_f(R
,rc
)*R_2
;
94 for(m
=0; (m
<DIM
); m
++) {
95 f_sr
[i
][m
] += dx
[m
]*fscal
;
96 f_sr
[j
][m
] -= dx
[m
]*fscal
;
103 fprintf(log
,"There were %d short range interactions, vsr=%g\n",ni
,vsr
);
108 void calc_ener(FILE *fp
,char *title
,gmx_bool bHeader
,int nmol
,
109 int natoms
,real phi
[],real charge
[],t_block
*excl
)
112 real qq
,qi
,vv
,V
,Vex
,Vc
,Vt
;
116 for(i
=0; (i
<natoms
); i
++) {
117 vv
= charge
[i
]*phi
[i
];
119 qq
+= charge
[i
]*charge
[i
];
122 Vc
= 0.5*C
*ONE_4PI_EPS0
*qq
;
125 for(i
=0; (i
<excl
->nr
); i
++) {
127 i2
= excl
->index
[i
+1];
129 for(j
=i1
; (j
<i2
); j
++) {
135 Vex
= qq
*0.5*C
*ONE_4PI_EPS0
;
140 fprintf(fp
,"%12s %12s %12s %12s %12s %12s\n",
141 "","Vphi","Vself","Vexcl","Vtot","1Mol");
144 fprintf(fp
,"%12s %12.5e %12.5e %12.5e %12.5e %12.5e\n",
145 title
,V
,Vc
,Vex
,Vt
,Vt
/nmol
);
148 void write_pqr(char *fn
,t_atoms
*atoms
,rvec x
[],real phi
[],real dx
)
153 fp
=gmx_fio_fopen(fn
,"w");
154 for(i
=0; (i
<atoms
->nr
); i
++) {
155 rnr
=atoms
->atom
[i
].resnr
;
156 fprintf(fp
,"%-6s%5d %-4.4s%3.3s %c%4d %8.3f%8.3f%8.3f%6.2f%6.2f\n",
157 "ATOM",(i
+1),*atoms
->atomname
[i
],*atoms
->resname
[rnr
],' ',
159 10*(dx
+x
[i
][XX
]),10*x
[i
][YY
],10*(x
[i
][ZZ
]),0.0,phi
[i
]);
164 void write_grid_pqr(char *fn
,int nx
,int ny
,int nz
,real
***phi
)
170 fp
=gmx_fio_fopen(fn
,"w");
171 for(i
=0; (i
<nx
); i
++)
172 for(j
=0; (j
<ny
); j
++)
173 for(k
=0; (k
<nz
); k
++,rnr
++)
174 fprintf(fp
,"%-6s%5d %-4.4s%3.3s %c%4d %8.3f%8.3f%8.3f%6.2f%6.2f\n",
175 "ATOM",(i
+1),"C","C",' ',
176 1+((rnr
+1) % 10000),fac
*i
,fac
*j
,fac
*k
,0.0,phi
[i
][j
][k
]);
179 void plot_phi(char *fn
,rvec box
,int natoms
,rvec x
[],real phi
[])
182 real phi_max
,rr
,gg
,bb
,fac
,dx
,x0
,y0
;
188 for(i
=0; (i
<natoms
); i
++)
189 phi_max
=max(phi_max
,fabs(phi
[i
]));
192 fprintf(stderr
,"All values zero, see .out file\n");
198 fprintf(stderr
,"Scaling box by %g\n",fac
);
201 (real
)(fac
*box
[XX
]+2*offset
),(real
)(fac
*box
[YY
]+2*offset
));
202 ps_translate(eps
,offset
,offset
);
204 ps_box(eps
,1,1,(real
)(fac
*box
[XX
]-1),(real
)(fac
*box
[YY
]-1));
206 for(i
=0; (i
<natoms
); i
++) {
209 gg
=bb
=(1.0+(phi
[i
]/phi_max
));
211 rr
=gg
=(1.0-(phi
[i
]/phi_max
));
215 ps_color(eps
,rr
,gg
,bb
);
218 ps_fillbox(eps
,(real
)(x0
-dx
),(real
)(y0
-dx
),(real
)(x0
+dx
),(real
)(y0
+dx
));
223 void plot_qtab(char *fn
,int nx
,int ny
,int nz
,real
***qtab
)
234 npt
=(box
[XX
]*box
[YY
]*box
[ZZ
]);
241 for(ix
=-nx
; (ix
<nx
); ix
++)
242 for(iy
=-ny
; (iy
<ny
); iy
++)
243 for(iz
=-nz
; (iz
<nz
); iz
++,i
++) {
246 xx
[i
][ZZ
]=iz
+nz
+0.5; /* onzin */
247 phi
[i
]=qtab
[ix
+nx
][iy
+ny
][iz
+nz
];
250 plot_phi(fn
,box
,npt
,xx
,phi
);
256 void print_phi(char *fn
,int natoms
,rvec x
[],real phi
[])
261 fp
=gmx_fio_fopen(fn
,"w");
262 for(i
=0; (i
<natoms
); i
++)
263 fprintf(fp
,"%10d %12.5e\n",i
,phi
[i
]);
267 void pr_f(char *fn
,int natoms
,rvec f
[])
272 fp
=gmx_fio_fopen(fn
,"w");
273 for(i
=0; (i
<natoms
); i
++)
274 fprintf(fp
," %12.5e\n %12.5e\n %12.5e\n",f
[i
][XX
],f
[i
][YY
],f
[i
][ZZ
]);
278 void test_pppm(FILE *log
, gmx_bool bVerbose
,
279 gmx_bool bGenerGhat
, char *ghatfn
,
280 t_atoms
*atoms
, t_inputrec
*ir
,
282 real charge
[], rvec box
,
283 real phi
[], real phi_s
[],
284 int nmol
, t_commrec
*cr
,
285 gmx_bool bOld
, t_block
*cgs
)
294 calc_nsb(cgs
,1,&nsb
,0);
296 /* First time only setup is done! */
297 init_pppm(log
,cr
,&nsb
,bVerbose
,bOld
,box
,ghatfn
,ir
);
299 ener
= do_pppm(log
,bVerbose
,x
,f
,charge
,box
,phi
,cr
,&nsb
,&nrnb
);
300 fprintf(log
,"Vpppm = %g\n",ener
);
302 sprintf(buf
,"PPPM-%d.pdb",ir
->nkx
);
303 write_pqr(buf
,atoms
,x
,phi
,0);
305 pr_f("pppm-force",atoms
->nr
,f
);
307 calc_ener(log
,buf
,FALSE
,nmol
,atoms
->nr
,phi
,charge
,&atoms
->excl
);
309 for(i
=0; (i
<atoms
->nr
); i
++)
311 sprintf(buf
,"PPPM-%d+SR",ir
->nkx
);
312 calc_ener(log
,buf
,FALSE
,nmol
,atoms
->nr
,phi
,charge
,&atoms
->excl
);
314 write_pqr(buf
,atoms
,x
,phi
,0);
317 void test_poisson(FILE *log
, gmx_bool bVerbose
,
318 t_atoms
*atoms
, t_inputrec
*ir
,
320 real charge
[], rvec box
,
321 real phi
[], real phi_s
[],
322 int nmol
, t_commrec
*cr
,
323 gmx_bool bFour
, rvec f_four
[],
324 real phi_f
[], gmx_bool bOld
)
334 /* First time only setup is done! */
336 for(i
=0; (i
<atoms
->nr
); i
++)
337 phi_f
[i
] -= phi_s
[i
];
338 ener
= do_optimize_poisson(log
,bVerbose
,ir
,atoms
->nr
,x
,f
,charge
,
339 box
,phi
,cr
,&nrnb
,f_four
,phi_f
,beta
,bOld
);
340 for(i
=0; (i
<atoms
->nr
); i
++)
341 phi_f
[i
] += phi_s
[i
];
345 ener
= do_poisson(log
,bVerbose
,ir
,atoms
->nr
,x
,f
,charge
,box
,phi
,
349 fprintf(log
,"Vpoisson = %g, nit = %d\n",ener
,nit
);
351 sprintf(buf
,"POISSON-%d.pdb",ir
->nkx
);
352 write_pqr(buf
,atoms
,x
,phi
,0);
354 pr_f("poisson-force",atoms
->nr
,f
);
356 calc_ener(log
,buf
,FALSE
,nmol
,atoms
->nr
,phi
,charge
,&atoms
->excl
);
358 for(i
=0; (i
<atoms
->nr
); i
++)
360 sprintf(buf
,"POISSON-%d+SR",ir
->nkx
);
361 calc_ener(log
,buf
,FALSE
,nmol
,atoms
->nr
,phi
,charge
,&atoms
->excl
);
363 write_pqr(buf
,atoms
,x
,phi
,0);
366 void test_four(FILE *log
,int NFILE
,t_filenm fnm
[],t_atoms
*atoms
,
367 t_inputrec
*ir
,rvec x
[],rvec f
[],rvec box
,real charge
[],
368 real phi_f
[],real phi_s
[],int nmol
,t_commrec
*cr
,
369 gmx_bool bOld
,gmx_bool bOldEwald
)
375 init_ewald_tab(&et
, NULL
, ir
, log
);
378 energy
= do_ewald(log
,ir
,atoms
->nr
,x
,f
,charge
,box
,phi_f
,cr
,bOld
,et
);
380 energy
= do_ewald_new(log
,ir
,atoms
->nr
,x
,f
,charge
,box
,phi_f
,cr
,bOld
,et
);
382 /*symmetrize_phi(log,atoms->nr,phi_f,bVerbose);*/
384 /* Plot the long range fourier solution in a matrix */
385 plot_phi(opt2fn("-elf",NFILE
,fnm
),box
,atoms
->nr
,x
,phi_f
);
386 print_phi(opt2fn("-of",NFILE
,fnm
),atoms
->nr
,x
,phi_f
);
387 calc_ener(log
,"Fourier",FALSE
,nmol
,atoms
->nr
,phi_f
,charge
,&atoms
->excl
);
388 write_pqr(opt2fn("-pf",NFILE
,fnm
),atoms
,x
,phi_f
,0.0/*1.5*box[XX]*/);
389 pr_f("four-force",atoms
->nr
,f
);
391 /* Calc and plot the total potential */
392 for(i
=0; (i
<atoms
->nr
); i
++) {
394 /*clear_rvec(f[i]);*/
396 calc_ener(log
,"Fourier+SR",FALSE
,nmol
,atoms
->nr
,phi_f
,charge
,&atoms
->excl
);
399 static void print_opts(FILE *fp
,t_inputrec
*ir
,gmx_bool bFour
)
401 fprintf(fp
,"Ewald solution: %s\n",gmx_bool_names
[bFour
]);
402 fprintf(fp
,"r1: %10.3e\n",ir
->rcoulomb_switch
);
403 fprintf(fp
,"rc: %10.3e\n",ir
->rcoulomb
);
405 fprintf(fp
,"KVectors: %10d %10d %10d\n",ir
->nkx
,ir
->nky
,ir
->nkz
);
409 int main(int argc
,char *argv
[])
411 static char *desc
[] = {
412 "[TT]testlr[tt] tests the PPPM and Ewald method for the",
413 "long range electrostatics problem."
415 static t_filenm fnm
[] = {
416 { efTPX
, NULL
, NULL
, ffREAD
},
417 { efHAT
, "-g", "ghat", ffOPTRD
},
418 { efOUT
, "-o", "rho", ffOPTWR
},
419 { efOUT
, "-op", "lr-pb", ffOPTWR
},
420 { efOUT
, "-of", "lr-four", ffOPTWR
},
421 { efOUT
, "-opt", "tot-pb", ffOPTWR
},
422 { efOUT
, "-oft", "tot-four", ffOPTWR
},
423 { efOUT
, "-fin", "lr-four", ffOPTWR
},
424 { efEPS
, "-es", "sr", ffOPTWR
},
425 { efEPS
, "-elf", "lr-four", ffOPTWR
},
426 { efEPS
, "-etf", "tot-four", ffOPTWR
},
427 { efEPS
, "-qr", "qk-real", ffOPTWR
},
428 { efEPS
, "-qi", "qk-im", ffOPTWR
},
429 { efEPS
, "-elp", "lr-pb", ffOPTWR
},
430 { efEPS
, "-etp", "tot-pb", ffOPTWR
},
431 { efEPS
, "-rho", "rho", ffOPTWR
},
432 { efEPS
, "-qq", "charge", ffOPTWR
},
433 { efXVG
, "-gt", "gk-tab", ffOPTWR
},
434 { efXVG
, "-fcorr","fcorr", ffWRITE
},
435 { efXVG
, "-pcorr","pcorr", ffWRITE
},
436 { efXVG
, "-ftotcorr","ftotcorr", ffWRITE
},
437 { efXVG
, "-ptotcorr","ptotcorr", ffWRITE
},
438 { efLOG
, "-l", "fptest", ffWRITE
},
439 { efXVG
, "-gr", "spread", ffOPTWR
},
440 { efPDB
, "-pf", "pqr-four", ffOPTWR
},
441 { efPDB
, "-phitot", "pppm-phitot", ffOPTWR
}
443 #define NFILE asize(fnm)
453 int i
,step
,nre
,natoms
,nmol
;
454 rvec
*x
,*f_sr
,*f_excl
,*f_four
,*f_pppm
,*f_pois
,box_size
,hbox
;
456 real t
,lambda
,vsr
,*charge
,*phi_f
,*phi_pois
,*phi_s
,*phi_p3m
,*rho
;
459 static gmx_bool bFour
=FALSE
,bVerbose
=FALSE
,bGGhat
=FALSE
,bPPPM
=TRUE
,
460 bPoisson
=FALSE
,bOld
=FALSE
,bOldEwald
=TRUE
;
461 static int nprocs
= 1;
462 static t_pargs pa
[] = {
463 { "-np", FALSE
, etINT
, &nprocs
, "Do it in parallel" },
464 { "-ewald", FALSE
, etBOOL
, &bFour
, "Do an Ewald solution"},
465 { "-pppm", FALSE
, etBOOL
, &bPPPM
, "Do a PPPM solution" },
466 { "-poisson",FALSE
, etBOOL
, &bPoisson
,"Do a Poisson solution" },
467 { "-v", FALSE
, etBOOL
, &bVerbose
,"Verbose on"},
468 { "-ghat", FALSE
, etBOOL
, &bGGhat
, "Generate Ghat function"},
469 { "-old", FALSE
, etBOOL
, &bOld
, "Use old function types"},
470 { "-oldewald",FALSE
,etBOOL
, &bOldEwald
,"Use old Ewald code"}
473 CopyRight(stderr
,argv
[0]);
474 parse_common_args_r(&argc
,argv
,PCA_CAN_TIME
| PCA_CAN_VIEW
,
475 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),desc
,0,NULL
,&oenv
);
478 cr
= init_par(&argc
,argv
);
479 open_log(ftp2fn(efLOG
,NFILE
,fnm
),cr
);
483 cr
= init_par(&argc
,argv
);
484 log
= ftp2FILE(efLOG
,NFILE
,fnm
,"w");
488 /* Read topology and coordinates */
489 read_tpxheader(ftp2fn(efTPX
,NFILE
,fnm
),&stath
,FALSE
);
490 snew(x
,stath
.natoms
);
491 snew(f_sr
,stath
.natoms
);
492 snew(f_excl
,stath
.natoms
);
493 snew(f_four
,stath
.natoms
);
494 snew(f_pppm
,stath
.natoms
);
495 snew(f_pois
,stath
.natoms
);
496 read_tpx(ftp2fn(efTPX
,NFILE
,fnm
),&step
,&t
,&lambda
,&ir
,
497 box
,&natoms
,x
,NULL
,NULL
,&top
);
498 excl
=&(top
.atoms
.excl
);
499 nmol
=top
.blocks
[ebMOLS
].nr
;
501 /* Allocate space for potential, charges and rho (charge density) */
502 snew(charge
,stath
.natoms
);
503 snew(phi_f
,stath
.natoms
);
504 snew(phi_p3m
,stath
.natoms
);
505 snew(phi_pois
,stath
.natoms
);
506 snew(phi_s
,stath
.natoms
);
507 snew(rho
,stath
.natoms
);
509 /* Set the charges */
510 for(i
=0; (i
<natoms
); i
++)
511 charge
[i
]=top
.atoms
.atom
[i
].q
;
513 /* Make a simple box vector instead of tensor */
514 for(i
=0; (i
<DIM
); i
++)
515 box_size
[i
]=box
[i
][i
];
517 /* Set some constants */
519 mdatoms
= atoms2md(&(top
.atoms
),FALSE
,FALSE
);
521 set_LRconsts(log
,ir
.rcoulomb_switch
,ir
.rcoulomb
,box_size
,fr
);
522 init_forcerec(log
,fr
,&ir
,&(top
.blocks
[ebMOLS
]),cr
,
523 &(top
.blocks
[ebCGS
]),&(top
.idef
),mdatoms
,box
,FALSE
);
524 calc_shifts(box
,box_size
,fr
->shift_vec
,FALSE
);
526 /* Periodicity stuff */
527 graph
= mk_graph(&(top
.idef
),top
.atoms
.nr
,FALSE
,FALSE
);
528 shift_self(graph
,fr
->shift_vec
,x
);
530 calc_LRcorrections(log
,0,natoms
,ir
.rcoulomb_switch
,
531 ir
.rcoulomb
,charge
,excl
,x
,f_excl
,bOld
);
532 pr_f("f_excl.dat",natoms
,f_excl
);
534 /* Compute the short range potential */
535 put_atoms_in_box(natoms
,box
,x
);
536 vsr
=phi_sr(log
,natoms
,x
,charge
,ir
.rcoulomb
,
537 ir
.rcoulomb_switch
,box_size
,phi_s
,excl
,f_sr
,bOld
);
538 pr_f("f_sr.dat",natoms
,f_sr
);
540 /* Plot the short range potential in a matrix */
541 calc_ener(log
,"Short Range",TRUE
,nmol
,natoms
,phi_s
,charge
,excl
);
545 test_four(log
,NFILE
,fnm
,&(top
.atoms
),&ir
,x
,f_four
,box_size
,charge
,phi_f
,
546 phi_s
,nmol
,cr
,bOld
,bOldEwald
);
549 test_pppm(log
,bVerbose
,bGGhat
,opt2fn("-g",NFILE
,fnm
),
550 &(top
.atoms
),&ir
,x
,f_pppm
,charge
,box_size
,phi_p3m
,phi_s
,nmol
,
551 cr
,bOld
,&(top
.blocks
[ebCGS
]));
554 test_poisson(log
,bVerbose
,
555 &(top
.atoms
),&ir
,x
,f_pois
,charge
,box_size
,phi_pois
,
556 phi_s
,nmol
,cr
,bFour
,f_four
,phi_f
,bOld
);
559 analyse_diff(log
,"PPPM",oenv
,
560 top
.atoms
.nr
,f_four
,f_pppm
,phi_f
,phi_p3m
,phi_s
,
561 opt2fn("-fcorr",NFILE
,fnm
),
562 opt2fn("-pcorr",NFILE
,fnm
),
563 opt2fn("-ftotcorr",NFILE
,fnm
),
564 opt2fn("-ptotcorr",NFILE
,fnm
));
566 if (bPoisson
&& bFour
)
567 analyse_diff(log
,"Poisson",oenv
,
568 top
.atoms
.nr
,f_four
,f_pois
,phi_f
,phi_pois
,phi_s
,
569 opt2fn("-fcorr",NFILE
,fnm
),
570 opt2fn("-pcorr",NFILE
,fnm
),
571 opt2fn("-ftotcorr",NFILE
,fnm
),
572 opt2fn("-ptotcorr",NFILE
,fnm
));