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
55 #include "thread_mpi.h"
57 #define p2(x) ((x)*(x))
58 #define p3(x) ((x)*(x)*(x))
59 #define p4(x) ((x)*(x)*(x)*(x))
61 static real A,A_3,B,B_4,C,c1,c2,c3,c4,c5,c6,One_4pi,FourPi_V,Vol,N0;
63 static tMPI_Thread_mutex_t shift_mutex=TMPI_THREAD_MUTEX_INITIALIZER;
67 void set_shift_consts(FILE *log,real r1,real rc,rvec box,t_forcerec *fr)
70 /* at the very least we shouldn't allow multiple threads to set these
72 tMPI_Thread_mutex_lock(&shift_mutex);
74 /* A, B and C are recalculated in tables.c */
76 A = (2*r1-5*rc)/(p3(rc)*p2(rc-r1));
77 B = (4*rc-2*r1)/(p3(rc)*p3(rc-r1));
78 /*C = (10*rc*rc-5*rc*r1+r1*r1)/(6*rc*rc); Hermans Eq. not correct */
81 gmx_fatal(FARGS,"r1 (%f) >= rc (%f) in %s, line %d",
82 r1,rc,__FILE__,__LINE__);
86 C = 1/rc-A_3*p3(rc-r1)-B_4*p4(rc-r1);
87 N0 = 2.0*M_PI*p3(rc)*p3(rc-r1);
89 Vol =(box[XX]*box[YY]*box[ZZ]);
90 FourPi_V=4.0*M_PI/Vol;
93 fprintf(debug,"Constants for short-range and fourier stuff:\n"
94 "r1 = %10.3f, rc = %10.3f\n"
95 "A = %10.3e, B = %10.3e, C = %10.3e, FourPi_V = %10.3e\n",
96 r1,rc,A,B,C,FourPi_V);
99 /* Constants derived by Mathematica */
100 c1 = -40*rc*rc + 50*rc*r1 - 16*r1*r1;
102 c3 = -10*rc*rc*rc + 20*rc*rc*r1 - 13*rc*r1*r1 + 3*r1*r1*r1;
103 c4 = -20*rc*rc + 40*rc*r1 - 14*r1*r1;
105 c6 = -5*rc*rc*r1 + 7*rc*r1*r1 - 2*r1*r1*r1;
108 fprintf(debug,"c1 = %10.3e, c2 = %10.3e, c3 = %10.3e\n"
109 "c4 = %10.3e, c5 = %10.3e, c6 = %10.3e, N0 = %10.3e\n",
110 c1,c2,c3,c4,c5,c6,N0);
113 One_4pi = 1.0/(4.0*M_PI);
115 tMPI_Thread_mutex_unlock(&shift_mutex);
119 real gk(real k,real rc,real r1)
120 /* Spread function in Fourier space. Eqs. 56-64 */
126 /* c1 thru c6 consts are global variables! */
127 gg = (FourPi_V / (N*k)) *
128 ( c1*k*cos(k*rc) + (c2+c3*k*k)*sin(k*rc) +
129 c4*k*cos(k*r1) + (c5+c6*k*k)*sin(k*r1) );
134 real gknew(real k,real rc,real r1)
141 return -15.0*((rck2-3.0)*sin(rck) + 3*rck*cos(rck))/(Vol*rck2*rck2*rck);
144 real calc_dx2dx(rvec xi,rvec xj,rvec box,rvec dx)
150 for(m=0; (m<DIM); m++) {
152 if (ddx < -0.5*box[m])
154 else if (ddx >= 0.5*box[m])
162 real calc_dx2(rvec xi,rvec xj,rvec box)
166 return calc_dx2dx(xi,xj,box,dx);
169 real shiftfunction(real r1,real rc,real R)
180 return A*dr*dr+B*dr*dr*dr;
183 real new_f(real r,real rc)
190 return 1.5*rrc2*rrc3 - 2.5*rrc3 + 1.0;
193 real new_phi(real r,real rc)
199 return 1/r-(0.125/rc)*(15 + 3*rrr*rrr - 10*rrr);
202 real old_f(real r,real rc,real r1)
213 return 1+A*r2*dr*dr+B*r2*dr*dr*dr;
216 real old_phi(real r,real rc,real r1)
227 return 1/r-A_3*dr*dr*dr-B_4*dr*dr*dr*dr - C;
230 real spreadfunction(real r1,real rc,real R)
241 return -One_4pi*(dr/R)*(2*A*(2*R-r1)+B*dr*(5*R-2*r1));
244 real potential(real r1,real rc,real R)
249 return (1.0/R - A_3*p3(R-r1) - B_4*p4(R-r1) - C);
256 real shift_LRcorrection(FILE *fp,int start,int natoms,
257 t_commrec *cr,t_forcerec *fr,
258 real charge[],t_blocka *excl,rvec x[],
259 gmx_bool bOld,matrix box,matrix lr_vir)
261 static gmx_bool bFirst=TRUE;
263 int i,i1,i2,j,k,m,iv,jv;
265 double qq; /* Necessary for precision */
266 double isp=0.564189583547756;
267 real qi,dr,dr2,dr_1,dr_3,fscal,Vexcl,qtot=0;
269 real r1=fr->rcoulomb_switch;
270 real rc=fr->rcoulomb;
275 for(i=start; (i<start+natoms); i++)
276 qq += charge[i]*charge[i];
278 /* Obsolete, until we implement dipole and charge corrections.
280 for(i=0;i<nsb->natoms;i++)
284 Vself = 0.5*C*ONE_4PI_EPS0*qq;
286 fprintf(fp,"Long Range corrections for shifted interactions:\n");
287 fprintf(fp,"r1 = %g, rc=%g\n",r1,rc);
288 fprintf(fp,"start=%d,natoms=%d\n",start,natoms);
289 fprintf(fp,"qq = %g, Vself=%g\n",qq,Vself);
296 for(i=start; (i<start+natoms); i++) {
297 /* Initiate local variables (for this i-particle) to 0 */
299 i2 = excl->index[i+1];
300 qi = charge[i]*ONE_4PI_EPS0;
302 /* Loop over excluded neighbours */
303 for(j=i1; (j<i2); j++) {
306 * First we must test whether k <> i, and then, because the
307 * exclusions are all listed twice i->k and k->i we must select
308 * just one of the two.
309 * As a minor optimization we only compute forces when the charges
316 rvec_sub(x[i],x[k],dx);
317 for(m=DIM-1; m>=0; m--) {
318 if (dx[m] > 0.5*box[m][m])
320 else if (dx[m] < -0.5*box[m][m])
325 dr_1 = gmx_invsqrt(dr2);
327 dr_3 = dr_1*dr_1*dr_1;
328 /* Compute exclusion energy and scalar force */
330 Vexcl += qq*(dr_1-potential(r1,rc,dr));
331 fscal = qq*(-shiftfunction(r1,rc,dr))*dr_3;
333 if ((fscal != 0) && debug)
334 fprintf(debug,"i: %d, k: %d, dr: %.3f fscal: %.3f\n",i,k,dr,fscal);
336 /* The force vector is obtained by multiplication with the
340 rvec_inc(fr->f_novirsum[k],df);
341 rvec_dec(fr->f_novirsum[i],df);
342 for(iv=0;iv<DIM;iv++)
343 for(jv=0;jv<DIM;jv++)
344 lr_vir[iv][jv]+=0.5*dx[iv]*df[jv];
350 fprintf(fp,"Long Range correction: Vexcl=%g\n",Vexcl);
353 /* Return the correction to the energy */
354 return (-(Vself+Vexcl));
357 real phi_aver(int natoms,real phi[])
363 for(i=0; (i<natoms); i++)
364 phitot=phitot+phi[i];
366 return (phitot/natoms);
369 real symmetrize_phi(FILE *log,int natoms,real phi[],gmx_bool bVerbose)
374 phitot=phi_aver(natoms,phi);
376 fprintf(log,"phi_aver = %10.3e\n",phitot);
378 for(i=0; (i<natoms); i++)
384 static real rgbset(real col)
394 real analyse_diff(FILE *log,char *label,const output_env_t oenv,
395 int natom,rvec ffour[],rvec fpppm[],
396 real phi_f[],real phi_p[],real phi_sr[],
397 char *fcorr,char *pcorr,char *ftotcorr,char *ptotcorr)
398 /* Analyse difference between forces from fourier (_f) and other (_p)
399 * LR solvers (and potential also).
400 * If the filenames are given, xvgr files are written.
401 * returns the root mean square error in the force.
405 FILE *fp=NULL,*gp=NULL;
406 real f2sum=0,p2sum=0;
407 real df,fmax,dp,pmax,rmsf;
409 fmax = fabs(ffour[0][0]-fpppm[0][0]);
410 pmax = fabs(phi_f[0] - phi_p[0]);
412 for(i=0; (i<natom); i++) {
413 dp = fabs(phi_f[i] - phi_p[i]);
416 for(m=0; (m<DIM); m++) {
417 df = fabs(ffour[i][m] - fpppm[i][m]);
423 rmsf = sqrt(f2sum/(3.0*natom));
424 fprintf(log,"\n********************************\nERROR ANALYSIS for %s\n",
426 fprintf(log,"%-10s%12s%12s\n","Error:","Max Abs","RMS");
427 fprintf(log,"%-10s %10.3f %10.3f\n","Force",fmax,rmsf);
428 fprintf(log,"%-10s %10.3f %10.3f\n","Potential",pmax,sqrt(p2sum/(natom)));
431 fp = xvgropen(fcorr,"LR Force Correlation","Four-Force","PPPM-Force",oenv);
432 for(i=0; (i<natom); i++) {
433 for(m=0; (m<DIM); m++) {
434 fprintf(fp,"%10.3f %10.3f\n",ffour[i][m],fpppm[i][m]);
438 do_view(oenv,fcorr,NULL);
441 fp = xvgropen(pcorr,"LR Potential Correlation","Four-Pot","PPPM-Pot",oenv);
443 gp = xvgropen(ptotcorr,"Total Potential Correlation","Four-Pot","PPPM-Pot",
445 for(i=0; (i<natom); i++) {
447 fprintf(fp,"%10.3f %10.3f\n",phi_f[i],phi_p[i]);
449 fprintf(gp,"%10.3f %10.3f\n",phi_f[i]+phi_sr[i],phi_p[i]+phi_sr[i]);
453 do_view(oenv,pcorr,NULL);
457 do_view(oenv,ptotcorr,NULL);