3 * This source code is part of
7 * GROningen MAchine for Chemical Simulations
9 * VERSION 3.3.99_development_20071104
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-2006, 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 Machine for Chemical Simulation
42 #include "gromacs/commandline/pargs.h"
44 #include "gromacs/utility/fatalerror.h"
45 #include "gromacs/fileio/xvgr.h"
47 #include "gromacs/fileio/pdbio.h"
49 #include "gromacs/utility/smalloc.h"
50 #include "gromacs/math/vec.h"
52 #include "gromacs/math/units.h"
55 #include "gromacs/fileio/trnio.h"
57 #include "gromacs/fileio/confio.h"
59 real
pot(real x
,real qq
,real c6
,real c12
)
61 return c12
*pow(x
,-12)-c6
*pow(x
,-6)+qq
*ONE_4PI_EPS0
/x
;
64 real
dpot(real x
,real qq
,real c6
,real c12
)
66 return -(12*c12
*pow(x
,-13)-6*c6
*pow(x
,-7)+qq
*ONE_4PI_EPS0
/sqr(x
));
69 int main(int argc
,char *argv
[])
71 static char *desc
[] = {
74 static real c6
=1.0e-3,c12
=1.0e-6,qi
=1,qj
=2,sig
=0.3,eps
=1,sigfac
=0.7;
76 { "-c6", FALSE
, etREAL
, {&c6
}, "c6" },
77 { "-c12", FALSE
, etREAL
, {&c12
}, "c12" },
78 { "-sig", FALSE
, etREAL
, {&sig
}, "sig" },
79 { "-eps", FALSE
, etREAL
, {&eps
}, "eps" },
80 { "-qi", FALSE
, etREAL
, {&qi
}, "qi" },
81 { "-qj", FALSE
, etREAL
, {&qj
}, "qj" },
82 { "-sigfac", FALSE
, etREAL
, {&sigfac
}, "Factor in front of sigma for starting the plot" }
85 { efXVG
, "-o", "potje", ffWRITE
}
87 #define NFILE asize(fnm)
91 real qq
,x
,oldx
,minimum
,mval
,dp
[2],pp
[2];
95 /* CopyRight(stdout,argv[0]);*/
96 parse_common_args(&argc
,argv
,PCA_CAN_VIEW
,
97 NFILE
,fnm
,asize(pa
),pa
,asize(desc
),
100 if (opt2parg_bSet("-sig",asize(pa
),pa
) ||
101 opt2parg_bSet("-eps",asize(pa
),pa
)) {
102 c6
= 4*eps
*pow(sig
,6);
103 c12
= 4*eps
*pow(sig
,12);
105 else if ((c6
!= 0) && (c12
!= 0)) {
106 sig
= pow(c12
/c6
,1.0/6.0);
112 printf("c6 = %12.5e, c12 = %12.5e\n",c6
,c12
);
113 printf("sigma = %12.5f, epsilon = %12.5f\n",sig
,eps
);
116 fp
= xvgropen(ftp2fn(efXVG
,NFILE
,fnm
),"Potential","r (nm)","E (kJ/mol)");
122 for(i
=0; (i
<100); i
++) {
123 x
= sigfac
*sig
+sig
*i
*0.02;
124 dp
[next
] = dpot(x
,qq
,c6
,c12
);
125 fprintf(fp
,"%10g %10g %10g\n",x
,pot(x
,qq
,c6
,c12
),
127 if ((i
> 0) && (dp
[cur
]*dp
[next
] < 0)) {
128 minimum
= oldx
+ dp
[cur
]*(x
-oldx
)/(dp
[cur
]-dp
[next
]);
129 mval
= pot(minimum
,qq
,c6
,c12
);
130 /*fprintf(stdout,"dp[cur] = %g, dp[next] = %g oldx = %g, dx = %g\n",
131 dp[cur],dp[next],oldx,x-oldx);*/
132 printf("Minimum at r = %g (nm). Value = %g (kJ/mol)\n",
141 do_view(ftp2fn(efXVG
,NFILE
,fnm
),NULL
);