chemshell_qchem / blob
2 #include <stdio.h>
3 #include <stdlib.h>
4 #include <math.h>
5 #include <unistd.h>
7 #include "objects.h"
8 #include "tcl.h"
9 #include "chemsh.h"
10 #include "dutil.h"
11 #include "dfragment.h"
12 #include "dmatrix.h"
14 declare_gradcheck(interp)
15 Tcl_Interp *interp;
16 {
17 int C_gradcheck_point();
18 int C_gradcheck_init();
19 Tcl_CreateCommand(interp, "gradcheck_init",C_gradcheck_init, (ClientData) "gradcheck_init",
20 (Tcl_CmdDeleteProc *) NULL);
21 Tcl_CreateCommand(interp, "gradcheck_point",C_gradcheck_point, (ClientData) "gradcheck_point",
22 (Tcl_CmdDeleteProc *) NULL);
23 return 0;
24 }
26 /* gradcheck control paramters */
28 static int n;
29 static double *x;
30 static double *x0, *fc, *fc1, f0, del;
31 static int onepoint=0;
32 static int icoord, ncoord, start, flag;
33 static Frag frag;
34 static Frag tempc;
35 static Matrix tempg, tempe;
37 static ObjList lc=NULL, ltc=NULL, lte=NULL, ltg=NULL;
39 int C_gradcheck_init(clientData, interp, argc, argv)
40 ClientData clientData;
41 Tcl_Interp *interp;
42 int argc;
43 char *argv[];
44 {
45 const char *s;
47 int i, j, np[MAXDIM], nd, iret;
48 char retbuf[5];
49 struct atom_struct *at;
51 del = 1.0e-5;
52 ncoord=99999;
53 start=0;
54 onepoint=0;
56 if(C_parsearg(interp,clientData,"formula del mincoord maxcoord coords tempc tempg tempe",argc,argv) != TCL_OK)return TCL_ERROR;
58 if(GetOptVar(s,"formula")){
59 if(!strcmp(s,"onepoint")){
60 onepoint=1;
61 } else if(!strcmp(s,"twopoint")){
62 onepoint=0;
63 } else {
64 pop_module_name(interp);
65 Tcl_SetResult(interp, "error in gradcheck_init - bad formula= argument ", TCL_VOLATILE);
66 return TCL_ERROR;
67 }
68 }
69 if(GetOptVar(s,"del"))if(Tcl_GetDouble(interp,s,&del) != TCL_OK){
70 pop_module_name(interp);
71 Tcl_SetResult(interp, "error in gradcheck_init - bad del= argument ", TCL_VOLATILE);
72 return TCL_ERROR;
73 }else{
74 printf("del set to %f\n",del);
75 }
77 if(GetOptVar(s,"maxcoord"))if(Tcl_GetInt(interp,s,&ncoord) != TCL_OK){
78 Tcl_SetResult(interp, "error in gradcheck_init - bad maxcoord= argument ", TCL_VOLATILE);
79 return TCL_ERROR;
80 }else{
81 printf("maxcoord set to %d\n",ncoord);
82 }
84 if(GetOptVar(s,"mincoord"))if(Tcl_GetInt(interp,s,&start) != TCL_OK){
85 Tcl_SetResult(interp, "error in gradcheck_init - bad mincoord= argument ", TCL_VOLATILE);
86 return TCL_ERROR;
87 }else{
88 printf("mincoord set to %d\n",start);
89 }
91 GetArg(s,"coords","gradcheck");
92 lc = get_objlist(s,"fragment",CHEMSH_OBJ_EXISTS,0);
94 if(!lc){
95 printf("%s: bad coords specification %s\n", (char *) clientData,s);
96 return TCL_ERROR;
97 }
99 frag = (Frag) lc->object->data;
101 /* scratch space */
102 GetArg(s,"tempc","gradcheck");
103 ltc = get_objlist(s,"fragment",CHEMSH_OBJ_CREATE,CHEMSH_OBJ_VOLATILE);
104 tempc = (Frag) ltc->object->data;
106 GetArg(s,"tempg","gradcheck");
107 ltg = get_objlist(s,"matrix",CHEMSH_OBJ_CREATE,CHEMSH_OBJ_VOLATILE);
108 tempg = (Matrix) ltg->object->data;
110 GetArg(s,"tempe","gradcheck");
111 lte = get_objlist(s,"matrix",CHEMSH_OBJ_CREATE,CHEMSH_OBJ_VOLATILE);
112 tempe = (Matrix) lte->object->data;
114 /* set the energy and gradient arrays to the correct size */
116 np[0]=3;
117 np[1]=frag->natom;
118 nd=2;
119 MATRIX_allocate_memory(tempg,nd,np,MATRIX_type_flag("real"));
121 np[0]=1;
122 np[1]=1;
123 nd=2;
124 MATRIX_allocate_memory(tempe,nd,np,MATRIX_type_flag("real"));
126 /* control information and workspace */
128 n = frag->natom*3;
129 x = (double *) ckalloc(sizeof(double) * n);
130 x0 = (double *) ckalloc(sizeof(double) * n);
131 fc = (double *) ckalloc(sizeof(double) * n);
132 if(! onepoint)fc1 = (double *) ckalloc(sizeof(double) * n);
133 else fc1=NULL;
135 if(ncoord > n)ncoord = n;
137 printf("gradcheck parameters : del %14.6e n %d formula %s\n",del,n,
138 ( onepoint ? "onepoint" : "twopoint" ) );
140 /* pack the reference coordinate array */
141 at=frag->atoms;
142 for(j=0,i=0;i<frag->natom;i++){
143 x0[j++]=at[i].pos.x[0];
144 x0[j++]=at[i].pos.x[1];
145 x0[j++]=at[i].pos.x[2];
146 }
148 /* make x a copy */
149 for(i=0;i<n;i++)x[i]=x0[i];
151 /* the inital structure is the first one to pass
152 to the energy and gradient routines */
153 copy_frag(tempc,frag,1);
155 printf("gradcheck initialisation done\n");
157 icoord=-1;
158 iret=2;
160 sprintf(retbuf,"%d",iret);
161 pop_module_name(interp);
162 Tcl_SetResult(interp, retbuf, TCL_VOLATILE);
163 return TCL_OK;
164 }
165 int C_gradcheck_point(clientData, interp, argc, argv)
166 ClientData clientData;
167 Tcl_Interp *interp;
168 int argc;
169 char *argv[];
170 {
171 int i, j, natom, iret, atom;
172 char retbuf[5], cart;
173 double sum;
174 struct atom_struct *at;
175 double gg;
176 double *g0;
178 if (argc != 1){
179 Tcl_AppendResult(interp, "No arguments required for ", argv[0], NULL);
180 return TCL_ERROR;
181 }
183 /* printf(">>>>> gradcheck_point\n"); */
185 at=tempc->atoms;
186 natom = tempc->natom;
187 g0 = tempg->d;
189 /* check the stuff loaded ok.. */
190 if(tempe->np[0] != 1 || tempe->np[1] != 1 || tempe->type != DENSEREAL){
191 printf("matrix size error for energy \n");
192 return TCL_ERROR;
193 }
195 if(icoord==-1){
196 printf(" Forces at reference geometry: \n");
198 /* check the stuff loaded ok.. */
199 if(tempg->np[0] != 3 || tempg->np[1] != natom || tempg->type != DENSEREAL){
200 printf("matrix size error for gradients \n");
201 return TCL_ERROR;
202 }
204 flag = 0;
205 f0=tempe->d[0];
206 sum = 0.0;
207 for(j=0,i=0;i<natom;i++){
208 sum = sum + tempg->d[i*3+0]*tempg->d[i*3+0] +
209 tempg->d[i*3+1]*tempg->d[i*3+1] +
210 tempg->d[i*3+2]*tempg->d[i*3+2];
211 printf("%8s %12.4e %12.4e %12.4e\n",
212 at[i].aname,
213 tempg->d[i*3+0],
214 tempg->d[i*3+1],
215 tempg->d[i*3+2]);
216 }
217 icoord += start;
218 }else if(icoord < ncoord){
220 if(onepoint){
221 fc[icoord]=tempe->d[0];
222 } else {
223 switch (flag) {
224 case 0:
225 fc[icoord]=tempe->d[0];
226 break;
227 default:
228 fc1[icoord]=tempe->d[0];
229 }
230 }
231 }
233 /*printf("icoord ncoord %d %d\n",icoord,ncoord);*/
235 if(icoord == ncoord-1 && (onepoint || !flag )){
236 printf("Results of gradient check:\n\n");
237 printf("coord FD Analytic error ratio \n");
238 for(i=start;i<ncoord;i++){
240 atom = i/3 +1;
241 switch(i-3*atom+3){
242 case 0:
243 cart='x';
244 break;
245 case 1:
246 cart='y';
247 break;
248 case 2:
249 cart='z';
250 break;
251 }
253 if(onepoint){
254 if( fabs(g0[i]) < 1.0e-10 ){
255 gg = (fc[i]-f0)/del;
256 printf("%4d %c %13.4e %13.4e %13.4e (g=0.0) \n",atom,cart,gg,g0[i],gg-g0[i]);
257 }else{
258 gg = (fc[i]-f0)/del;
259 /* printf("check %12.4e %12.4e %12.4e\n",fabs(g0[i]),1.0e-10,gg/g0[i]);*/
260 printf("%4d %c %13.4e %13.4e %13.4e %9.5f \n ",atom,cart,gg,g0[i],gg-g0[i],gg/g0[i]);
261 }
262 }else{ /* two point formula */
263 if( fabs(g0[i]) < 1.0e-10 ){
264 gg = (fc1[i]-fc[i])/(2.0*del);
265 printf("%4d %c %13.4e %13.4e %13.4e (g=0.0) \n",atom,cart,gg,g0[i],gg-g0[i]);
266 }else{
267 gg = (fc1[i]-fc[i])/(2.0*del);
268 /* printf("check %12.4e %12.4e %12.4e\n",fabs(g0[i]),1.0e-10,gg/g0[i]);*/
269 printf("%4d %c %13.4e %13.4e %13.4e %9.5f \n",atom,cart,gg,g0[i],gg-g0[i],gg/g0[i]);
270 }
271 }
272 }
274 if(lc)rel_objlist(lc);
275 if(ltc)rel_objlist(ltc);
276 if(lte)rel_objlist(lte);
277 if(ltg)rel_objlist(ltg);
279 ckfree((char *) x);
280 ckfree((char *) x0);
281 ckfree((char *) fc);
282 if(fc1)ckfree((char *) fc1);
284 iret = 0;
285 }else{
286 if(onepoint){
287 icoord++;
288 for(i=0;i<n;i++)x[i]=x0[i];
289 x[icoord]+=del;
290 printf(" stepping coord %d\n",icoord);
291 }else { /* 2 point formula */
292 switch (flag){
293 case 0:
294 icoord ++;
295 for(i=0;i<n;i++)x[i]=x0[i];
296 x[icoord]+=del;
297 printf(" stepping coord + %d\n",icoord);
298 break;
299 default:
300 for(i=0;i<n;i++)x[i]=x0[i];
301 x[icoord]-=del;
302 printf(" stepping coord - %d\n",icoord);
303 break;
304 }
305 flag = !flag;
306 }
307 /* update the structure */
308 sum=0.0;
309 for(j=0,i=0;i<tempc->natom;i++){
310 sum = sum + (at[i].pos.x[0] - x[j])*(at[i].pos.x[0] - x[j]);
311 at[i].pos.x[0]=x[j++];
312 sum = sum + (at[i].pos.x[1] - x[j])*(at[i].pos.x[1] - x[j]);
313 at[i].pos.x[1]=x[j++];
314 sum = sum + (at[i].pos.x[2] - x[j])*(at[i].pos.x[2] - x[j]);
315 at[i].pos.x[2]=x[j++];
316 }
317 iret = 1;
318 }
320 if(iret < 0){
321 Tcl_SetResult(interp, "error in gradcheck_point", TCL_VOLATILE);
322 return TCL_ERROR;
323 }
324 sprintf(retbuf,"%d",iret);
325 Tcl_SetResult(interp, retbuf, TCL_VOLATILE);
326 return TCL_OK;
327 }
