1 | |
2 | |
3 | |
4 | |
5 | |
6 | |
7 | |
8 | |
9 | |
10 | |
11 | |
12 | |
13 | |
14 | |
15 | |
16 | |
17 | |
18 | |
19 | |
20 | |
21 | |
22 | |
23 | |
24 | |
25 | |
26 | |
27 | |
28 | |
29 | |
30 | |
31 | |
32 | |
33 | |
34 | |
35 | |
36 | |
37 | #ifdef HAVE_CONFIG_H1 |
38 | #include <config.h> |
39 | #endif |
40 | |
41 | #include <math.h> |
42 | #include <stdio.h> |
43 | |
44 | #include "typedefs.h" |
45 | #include "gromacs/commandline/pargs.h" |
46 | #include "gromacs/fileio/xvgr.h" |
47 | #include "viewit.h" |
48 | #include "gromacs/fileio/pdbio.h" |
49 | #include "macros.h" |
50 | #include "gromacs/utility/smalloc.h" |
51 | #include "gromacs/math/vec.h" |
52 | #include "pbc.h" |
53 | #include "physics.h" |
54 | #include "names.h" |
55 | #include "txtdump.h" |
56 | #include "gromacs/fileio/trnio.h" |
57 | #include "symtab.h" |
58 | #include "gromacs/fileio/confio.h" |
59 | |
60 | real pot(real x, real qq, real c6, real cn, int npow) |
61 | { |
62 | return cn*pow(x, -npow)-c6*pow(x, -6)+qq*ONE_4PI_EPS0((332.0636930*(4.184))*0.1)/x; |
63 | } |
64 | |
65 | real bhpot(real x, real A, real B, real C) |
66 | { |
67 | return A*exp(-B*x) - C*pow(x, -6.0); |
68 | } |
69 | |
70 | real dpot(real x, real qq, real c6, real cn, int npow) |
71 | { |
72 | return -(npow*cn*pow(x, -npow-1)-6*c6*pow(x, -7)+qq*ONE_4PI_EPS0((332.0636930*(4.184))*0.1)/sqr(x)); |
73 | } |
74 | |
75 | int gmx_sigeps(int argc, char *argv[]) |
76 | { |
77 | const char *desc[] = { |
78 | "[THISMODULE] is a simple utility that converts C6/C12 or C6/Cn combinations", |
79 | "to [GRK]sigma[grk] and [GRK]epsilon[grk], or vice versa. It can also plot the potential", |
80 | "in file. In addition, it makes an approximation of a Buckingham potential", |
81 | "to a Lennard-Jones potential." |
82 | }; |
83 | static real c6 = 1.0e-3, cn = 1.0e-6, qi = 0, qj = 0, sig = 0.3, eps = 1, sigfac = 0.7; |
84 | static real Abh = 1e5, Bbh = 32, Cbh = 1e-3; |
85 | static int npow = 12; |
86 | t_pargs pa[] = { |
87 | { "-c6", FALSE0, etREAL, {&c6}, "C6" }, |
88 | { "-cn", FALSE0, etREAL, {&cn}, "Constant for repulsion" }, |
89 | { "-pow", FALSE0, etINT, {&npow}, "Power of the repulsion term" }, |
90 | { "-sig", FALSE0, etREAL, {&sig}, "[GRK]sigma[grk]" }, |
91 | { "-eps", FALSE0, etREAL, {&eps}, "[GRK]epsilon[grk]" }, |
92 | { "-A", FALSE0, etREAL, {&Abh}, "Buckingham A" }, |
93 | { "-B", FALSE0, etREAL, {&Bbh}, "Buckingham B" }, |
94 | { "-C", FALSE0, etREAL, {&Cbh}, "Buckingham C" }, |
95 | { "-qi", FALSE0, etREAL, {&qi}, "qi" }, |
96 | { "-qj", FALSE0, etREAL, {&qj}, "qj" }, |
97 | { "-sigfac", FALSE0, etREAL, {&sigfac}, "Factor in front of [GRK]sigma[grk] for starting the plot" } |
98 | }; |
99 | t_filenm fnm[] = { |
100 | { efXVG, "-o", "potje", ffWRITE1<<2 } |
101 | }; |
102 | output_env_t oenv; |
103 | #define NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))) asize(fnm)((int)(sizeof(fnm)/sizeof((fnm)[0]))) |
104 | const char *legend[] = { "Lennard-Jones", "Buckingham" }; |
105 | FILE *fp; |
106 | int i; |
107 | gmx_bool bBham; |
108 | real qq, x, oldx, minimum, mval, dp[2], pp[2]; |
109 | int cur = 0; |
110 | #define next(1-cur) (1-cur) |
111 | |
112 | if (!parse_common_args(&argc, argv, PCA_CAN_VIEW(1<<5), |
113 | NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm, asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa, asize(desc)((int)(sizeof(desc)/sizeof((desc)[0]))), |
114 | desc, 0, NULL((void*)0), &oenv)) |
115 | { |
116 | return 0; |
117 | } |
118 | |
119 | bBham = (opt2parg_bSet("-A", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa) || |
120 | opt2parg_bSet("-B", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa) || |
121 | opt2parg_bSet("-C", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa)); |
122 | |
123 | if (bBham) |
124 | { |
125 | c6 = Cbh; |
126 | sig = pow((6.0/npow)*pow(npow/Bbh, npow-6.0), 1.0/(npow-6.0)); |
127 | eps = c6/(4*pow(sig, 6.0)); |
128 | cn = 4*eps*pow(sig, npow); |
129 | } |
130 | else |
131 | { |
132 | if (opt2parg_bSet("-sig", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa) || |
133 | opt2parg_bSet("-eps", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa)) |
134 | { |
135 | c6 = 4*eps*pow(sig, 6); |
136 | cn = 4*eps*pow(sig, npow); |
137 | } |
138 | else if (opt2parg_bSet("-c6", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa) || |
139 | opt2parg_bSet("-cn", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa) || |
140 | opt2parg_bSet("-pow", asize(pa)((int)(sizeof(pa)/sizeof((pa)[0]))), pa)) |
141 | { |
142 | sig = pow(cn/c6, 1.0/(npow-6.0)); |
143 | eps = 0.25*c6*pow(sig, -6.0); |
144 | } |
145 | else |
146 | { |
147 | sig = eps = 0; |
148 | } |
149 | printf("c6 = %12.5e, c%d = %12.5e\n", c6, npow, cn); |
150 | printf("sigma = %12.5f, epsilon = %12.5f\n", sig, eps); |
151 | |
152 | minimum = pow(npow/6.0*pow(sig, npow-6.0), 1.0/(npow-6)); |
153 | printf("Van der Waals minimum at %g, V = %g\n\n", |
154 | minimum, pot(minimum, 0, c6, cn, npow)); |
155 | printf("Fit of Lennard Jones (%d-6) to Buckingham:\n", npow); |
156 | Bbh = npow/minimum; |
157 | Cbh = c6; |
158 | Abh = 4*eps*pow(sig/minimum, npow)*exp(npow); |
159 | printf("A = %g, B = %g, C = %g\n", Abh, Bbh, Cbh); |
160 | } |
161 | qq = qi*qj; |
162 | |
163 | fp = xvgropen(ftp2fn(efXVG, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), "Potential", "r (nm)", "E (kJ/mol)", |
164 | oenv); |
165 | xvgr_legend(fp, asize(legend)((int)(sizeof(legend)/sizeof((legend)[0]))), legend, |
166 | oenv); |
167 | if (sig == 0) |
168 | { |
169 | sig = 0.25; |
170 | } |
171 | minimum = -1; |
| Value stored to 'minimum' is never read |
172 | mval = 0; |
173 | oldx = 0; |
174 | for (i = 0; (i < 100); i++) |
175 | { |
176 | x = sigfac*sig+sig*i*0.02; |
177 | dp[next(1-cur)] = dpot(x, qq, c6, cn, npow); |
178 | fprintf(fp, "%10g %10g %10g\n", x, pot(x, qq, c6, cn, npow), |
179 | bhpot(x, Abh, Bbh, Cbh)); |
180 | if (qq != 0) |
181 | { |
182 | if ((i > 0) && (dp[cur]*dp[next(1-cur)] < 0)) |
183 | { |
184 | minimum = oldx + dp[cur]*(x-oldx)/(dp[cur]-dp[next(1-cur)]); |
185 | mval = pot(minimum, qq, c6, cn, npow); |
186 | printf("Van der Waals + Coulomb minimum at r = %g (nm). Value = %g (kJ/mol)\n", |
187 | minimum, mval); |
188 | } |
189 | } |
190 | cur = next(1-cur); |
191 | oldx = x; |
192 | |
193 | } |
194 | gmx_ffclose(fp); |
195 | |
196 | do_view(oenv, ftp2fn(efXVG, NFILE((int)(sizeof(fnm)/sizeof((fnm)[0]))), fnm), NULL((void*)0)); |
197 | |
198 | return 0; |
199 | } |