4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
12 * Copyright (c) 1991-2001
13 * BIOSON Research Institute, Dept. of Biophysical Chemistry
14 * University of Groningen, The Netherlands
16 * This program is free software; you can redistribute it and/or
17 * modify it under the terms of the GNU General Public License
18 * as published by the Free Software Foundation; either version 2
19 * of the License, or (at your option) any later version.
21 * If you want to redistribute modifications, please consider that
22 * scientific software is very special. Version control is crucial -
23 * bugs must be traceable. We will be happy to consider code for
24 * inclusion in the official distribution, but derived work must not
25 * be called official GROMACS. Details are found in the README & COPYING
26 * files - if they are missing, get the official version at www.gromacs.org.
28 * To help us fund GROMACS development, we humbly ask that you cite
29 * the papers on the package - you can find them in the top README file.
31 * Do check out http://www.gromacs.org , or mail us at gromacs@gromacs.org .
34 * Gromacs Runs On Most of All Computer Systems
36 static char *SRCID_sigeps_c = "$Id$";
56 real pot(real x,real qq,real c6,real c12)
58 return c12*pow(x,-12)-c6*pow(x,-6)+qq*ONE_4PI_EPS0/x;
61 real dpot(real x,real qq,real c6,real c12)
63 return -(12*c12*pow(x,-13)-6*c6*pow(x,-7)+qq*ONE_4PI_EPS0/sqr(x));
66 int main(int argc,char *argv[])
68 static char *desc[] = {
71 static real c6=1.0e-3,c12=1.0e-6,qi=1,qj=2,sig=0.3,eps=1,sigfac=0.7;
73 { "-c6", FALSE, etREAL, {&c6}, "c6" },
74 { "-c12", FALSE, etREAL, {&c12}, "c12" },
75 { "-sig", FALSE, etREAL, {&sig}, "sig" },
76 { "-eps", FALSE, etREAL, {&eps}, "eps" },
77 { "-qi", FALSE, etREAL, {&qi}, "qi" },
78 { "-qj", FALSE, etREAL, {&qj}, "qj" },
79 { "-sigfac", FALSE, etREAL, {&sigfac}, "Factor in front of sigma for starting the plot" }
82 { efXVG, "-o", "potje", ffWRITE }
84 #define NFILE asize(fnm)
88 real qq,x,oldx,minimum,mval,dp[2],pp[2];
92 /* CopyRight(stdout,argv[0]);*/
93 parse_common_args(&argc,argv,PCA_CAN_VIEW,
94 FALSE,NFILE,fnm,asize(pa),pa,asize(desc),
97 if (opt2parg_bSet("-sig",asize(pa),pa) ||
98 opt2parg_bSet("-eps",asize(pa),pa)) {
99 c6 = 4*eps*pow(sig,6);
100 c12 = 4*eps*pow(sig,12);
102 else if ((c6 != 0) && (c12 != 0)) {
103 sig = pow(c12/c6,1.0/6.0);
109 printf("c6 = %12.5e, c12 = %12.5e\n",c6,c12);
110 printf("sigma = %12.5f, epsilon = %12.5f\n",sig,eps);
113 fp = xvgropen(ftp2fn(efXVG,NFILE,fnm),"Potential","r (nm)","E (kJ/mol)");
119 for(i=0; (i<100); i++) {
120 x = sigfac*sig+sig*i*0.02;
121 dp[next] = dpot(x,qq,c6,c12);
122 fprintf(fp,"%10g %10g %10g\n",x,pot(x,qq,c6,c12),
124 if ((i > 0) && (dp[cur]*dp[next] < 0)) {
125 minimum = oldx + dp[cur]*(x-oldx)/(dp[cur]-dp[next]);
126 mval = pot(minimum,qq,c6,c12);
127 /*fprintf(stdout,"dp[cur] = %g, dp[next] = %g oldx = %g, dx = %g\n",
128 dp[cur],dp[next],oldx,x-oldx);*/
129 printf("Minimum at r = %g (nm). Value = %g (kJ/mol)\n",
138 do_view(ftp2fn(efXVG,NFILE,fnm),NULL);