4 * This source code is part of
8 * GROningen MAchine for Chemical Simulations
11 * Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
12 * Copyright (c) 1991-2000, University of Groningen, The Netherlands.
13 * Copyright (c) 2001-2004, The GROMACS development team,
14 * check out http://www.gromacs.org for more information.
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 * For more info, check our website at http://www.gromacs.org
34 * Good gRace! Old Maple Actually Chews Slate
55 real
pot(real x
,real qq
,real c6
,real cn
,int npow
)
57 return cn
*pow(x
,-npow
)-c6
*pow(x
,-6)+qq
*ONE_4PI_EPS0
/x
;
60 real
bhpot(real x
,real qq
,real A
,real B
,real C
)
62 return A
*exp(-B
*x
) - C
*pow(x
,-6.0);
65 real
dpot(real x
,real qq
,real c6
,real cn
,int npow
)
67 return -(npow
*cn
*pow(x
,-npow
-1)-6*c6
*pow(x
,-7)+qq
*ONE_4PI_EPS0
/sqr(x
));
70 int main(int argc
,char *argv
[])
72 static char *desc
[] = {
75 static real c6
=1.0e-3,cn
=1.0e-6,qi
=0,qj
=0,sig
=0.3,eps
=1,sigfac
=0.7;
76 static real Abh
=1e5
,Bbh
=32,Cbh
=1e-3;
79 { "-c6", FALSE
, etREAL
, {&c6
}, "c6" },
80 { "-cn", FALSE
, etREAL
, {&cn
}, "constant for repulsion" },
81 { "-pow", FALSE
, etINT
, {&npow
},"power of the repulsion term" },
82 { "-sig", FALSE
, etREAL
, {&sig
}, "sig" },
83 { "-eps", FALSE
, etREAL
, {&eps
}, "eps" },
84 { "-A", FALSE
, etREAL
, {&Abh
}, "Buckingham A" },
85 { "-B", FALSE
, etREAL
, {&Bbh
}, "Buckingham B" },
86 { "-C", FALSE
, etREAL
, {&Cbh
}, "Buckingham C" },
87 { "-qi", FALSE
, etREAL
, {&qi
}, "qi" },
88 { "-qj", FALSE
, etREAL
, {&qj
}, "qj" },
89 { "-sigfac", FALSE
, etREAL
, {&sigfac
}, "Factor in front of sigma for starting the plot" }
92 { efXVG
, "-o", "potje", ffWRITE
}
94 #define NFILE asize(fnm)
95 static char *legend
[] = { "Lennard-Jones", "Buckingham" };
99 real qq
,x
,oldx
,minimum
,mval
,dp
[2],pp
[2];
103 /* CopyRight(stdout,argv[0]);*/
104 parse_common_args(&argc
,argv
,PCA_CAN_VIEW
,
105 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),
108 bBham
= (opt2parg_bSet("-A",asize(pa
),pa
) ||
109 opt2parg_bSet("-B",asize(pa
),pa
) ||
110 opt2parg_bSet("-C",asize(pa
),pa
));
114 sig
= pow((6.0/npow
)*pow(npow
/Bbh
,npow
-6.0),1.0/(npow
-6.0));
115 eps
= c6
/(4*pow(sig
,6.0));
116 cn
= 4*eps
*pow(sig
,npow
);
119 if (opt2parg_bSet("-sig",asize(pa
),pa
) ||
120 opt2parg_bSet("-eps",asize(pa
),pa
)) {
121 c6
= 4*eps
*pow(sig
,6);
122 cn
= 4*eps
*pow(sig
,npow
);
124 else if (opt2parg_bSet("-c6",asize(pa
),pa
) ||
125 opt2parg_bSet("-cn",asize(pa
),pa
) ||
126 opt2parg_bSet("-pow",asize(pa
),pa
)) {
127 sig
= pow(cn
/c6
,1.0/(npow
-6.0));
128 eps
= 0.25*c6
*pow(sig
,-6.0);
133 printf("c6 = %12.5e, c%d = %12.5e\n",c6
,npow
,cn
);
134 printf("sigma = %12.5f, epsilon = %12.5f\n",sig
,eps
);
136 minimum
= pow(npow
/6.0*pow(sig
,npow
-6.0),1.0/(npow
-6));
137 printf("Van der Waals minimum at %g, V = %g\n\n",
138 minimum
,pot(minimum
,0,c6
,cn
,npow
));
139 printf("Fit of Lennard Jones (%d-6) to Buckingham:\n",npow
);
142 Abh
= 4*eps
*pow(sig
/minimum
,npow
)*exp(npow
);
143 printf("A = %g, B = %g, C = %g\n",Abh
,Bbh
,Cbh
);
147 fp
= xvgropen(ftp2fn(efXVG
,NFILE
,fnm
),"Potential","r (nm)","E (kJ/mol)");
148 xvgr_legend(fp
,asize(legend
),legend
);
154 for(i
=0; (i
<100); i
++) {
155 x
= sigfac
*sig
+sig
*i
*0.02;
156 dp
[next
] = dpot(x
,qq
,c6
,cn
,npow
);
157 fprintf(fp
,"%10g %10g %10g\n",x
,pot(x
,qq
,c6
,cn
,npow
),
158 bhpot(x
,qq
,Abh
,Bbh
,Cbh
));
160 if ((i
> 0) && (dp
[cur
]*dp
[next
] < 0)) {
161 minimum
= oldx
+ dp
[cur
]*(x
-oldx
)/(dp
[cur
]-dp
[next
]);
162 mval
= pot(minimum
,qq
,c6
,cn
,npow
);
163 printf("Van der Waals + Coulomb minimum at r = %g (nm). Value = %g (kJ/mol)\n",
173 do_view(ftp2fn(efXVG
,NFILE
,fnm
),NULL
);