4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Copyright (c) 1991-2001, University of Groningen, The Netherlands
12 * This program is free software; you can redistribute it and/or
13 * modify it under the terms of the GNU General Public License
14 * as published by the Free Software Foundation; either version 2
15 * of the License, or (at your option) any later version.
17 * If you want to redistribute modifications, please consider that
18 * scientific software is very special. Version control is crucial -
19 * bugs must be traceable. We will be happy to consider code for
20 * inclusion in the official distribution, but derived work must not
21 * be called official GROMACS. Details are found in the README & COPYING
22 * files - if they are missing, get the official version at www.gromacs.org.
24 * To help us fund GROMACS development, we humbly ask that you cite
25 * the papers on the package - you can find them in the top README file.
27 * For more info, check our website at http://www.gromacs.org
30 * Great Red Owns Many ACres of Sand
52 real
pot(real x
,real qq
,real c6
,real c12
)
54 return c12
*pow(x
,-12)-c6
*pow(x
,-6)+qq
*ONE_4PI_EPS0
/x
;
57 real
dpot(real x
,real qq
,real c6
,real c12
)
59 return -(12*c12
*pow(x
,-13)-6*c6
*pow(x
,-7)+qq
*ONE_4PI_EPS0
/sqr(x
));
62 int main(int argc
,char *argv
[])
64 static char *desc
[] = {
67 static real c6
=1.0e-3,c12
=1.0e-6,qi
=1,qj
=2,sig
=0.3,eps
=1,sigfac
=0.7;
69 { "-c6", FALSE
, etREAL
, {&c6
}, "c6" },
70 { "-c12", FALSE
, etREAL
, {&c12
}, "c12" },
71 { "-sig", FALSE
, etREAL
, {&sig
}, "sig" },
72 { "-eps", FALSE
, etREAL
, {&eps
}, "eps" },
73 { "-qi", FALSE
, etREAL
, {&qi
}, "qi" },
74 { "-qj", FALSE
, etREAL
, {&qj
}, "qj" },
75 { "-sigfac", FALSE
, etREAL
, {&sigfac
}, "Factor in front of sigma for starting the plot" }
78 { efXVG
, "-o", "potje", ffWRITE
}
80 #define NFILE asize(fnm)
84 real qq
,x
,oldx
,minimum
,mval
,dp
[2],pp
[2];
88 /* CopyRight(stdout,argv[0]);*/
89 parse_common_args(&argc
,argv
,PCA_CAN_VIEW
,
90 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),
93 if (opt2parg_bSet("-sig",asize(pa
),pa
) ||
94 opt2parg_bSet("-eps",asize(pa
),pa
)) {
95 c6
= 4*eps
*pow(sig
,6);
96 c12
= 4*eps
*pow(sig
,12);
98 else if ((c6
!= 0) && (c12
!= 0)) {
99 sig
= pow(c12
/c6
,1.0/6.0);
105 printf("c6 = %12.5e, c12 = %12.5e\n",c6
,c12
);
106 printf("sigma = %12.5f, epsilon = %12.5f\n",sig
,eps
);
109 fp
= xvgropen(ftp2fn(efXVG
,NFILE
,fnm
),"Potential","r (nm)","E (kJ/mol)");
115 for(i
=0; (i
<100); i
++) {
116 x
= sigfac
*sig
+sig
*i
*0.02;
117 dp
[next
] = dpot(x
,qq
,c6
,c12
);
118 fprintf(fp
,"%10g %10g %10g\n",x
,pot(x
,qq
,c6
,c12
),
120 if ((i
> 0) && (dp
[cur
]*dp
[next
] < 0)) {
121 minimum
= oldx
+ dp
[cur
]*(x
-oldx
)/(dp
[cur
]-dp
[next
]);
122 mval
= pot(minimum
,qq
,c6
,c12
);
123 /*fprintf(stdout,"dp[cur] = %g, dp[next] = %g oldx = %g, dx = %g\n",
124 dp[cur],dp[next],oldx,x-oldx);*/
125 printf("Minimum at r = %g (nm). Value = %g (kJ/mol)\n",
134 do_view(ftp2fn(efXVG
,NFILE
,fnm
),NULL
);