chemshell_qchem / blob
1 /******************************************************************************
3 COPYRIGHT - Copyright and any other intellectual property rights in this
4 software belong to CLRC. the recipient may not copy the software except for
5 backup purposes. Nor may the recipient re-distribute the software to any
6 third party without the prior written consent of CLRC.
8 WARRANTIES - This software is supplied on an "AS IS" basis and no
9 representations are made, nor warranties are given as to its accuracy,
10 usefulness or performance. Nor do CLRC give any warranties, express or
11 implied, as to its merchantability or fitness for a particular purpose.
13 LIABILITY - CLRC will not be liable for an loss of any kind arising from the
14 supply and use of this software.
16 ****************************************************************************
18 $Author: tbenighaus $
19 $Date: 2010-09-23 09:35:53 +0200 (Thu, 23. Sep 2010) $
20 $Revision: 3689 $
21 $Source$
22 $State$
24 ****************************************************************************
26 C code for the restraints function
28 Providing a flexible way to augment any energy/gradient
29 evaluator with some additional terms, typically for umbrella
30 sampling and related methods
32 */
34 #include <stdio.h>
35 #include <stdlib.h>
36 #include <unistd.h>
37 #include <sys/types.h>
38 #include <string.h>
39 #include <math.h>
41 #include "tcl.h"
42 #include "objects.h"
43 #include "chemsh.h"
44 #include "dutil.h"
45 #include "dblock.h"
46 #include "dmatrix.h"
47 #include "dfragment.h"
48 #include "restraint.h"
50 #if defined (_WIN32) || defined (CRAYXX)
51 #define CAPITALSYMBOLS
52 #endif
54 #ifdef CAPITALSYMBOLS
55 #define torsionval_ TORSIONVAL
56 #define torsionderivative_ TORSIONDERIVATIVE
57 #endif
59 #ifdef NOC_
60 #define torsionval_ torsionval
61 #define torsionderivative_ torsionderivative
62 #endif
64 extern Tcl_Interp *chemsh_interp;
66 static int debug = 0;
68 /* Declarations */
69 int rst_newconfig(Tcl_Interp *interp,struct rst_struct *rst, int argc,const char **argv);
70 int calc_rst_ff(Ff ff, Frag f, Matrix e, Matrix g);
71 int rst_set_config(Tcl_Interp *interp, struct rst_struct *rst);
72 Tcl_CmdProc Chemsh_RestraintObjCmd, RestraintCmd, DispersionCmd;
73 void torsionval_(double torpos[4][3], double *phi);
74 void torsionderivative_(double torpos[4][3], double *phi, double torder[4][3]);
77 void
78 declare_restraint(Tcl_Interp *interp){
79 Tcl_CreateCommand(interp, "restraint_c" ,RestraintCmd,
80 (ClientData) NULL,(Tcl_CmdDeleteProc *) NULL);
81 Tcl_CreateCommand(interp, "dispersion_c" ,DispersionCmd,
82 (ClientData) NULL,(Tcl_CmdDeleteProc *) NULL);
83 }
86 int
87 RestraintCmd(ClientData clientData,
88 Tcl_Interp *interp,
89 int argc,
90 const char *argv[])
91 {
93 int Chemsh_RestraintObjCmd();
95 /* Local variables */
96 struct rst_struct *rst;
98 const char *s, **work, **work2;
99 int itmp, i, j;
100 double atmp;
101 const char *cmd;
102 int nfld;
103 Ff ff;
104 int ff_counter;
106 /* string-buffers for use in Tcl_Eval , Tcl_SetVar
107 (to avoid conflict with strbuff altered inside the called C-routines */
108 char tclbuff[STRBUFFLEN];
109 char tclbuff2[STRBUFFLEN];
111 char buff2[2];
113 cmd = argv[0];
114 rst = (struct rst_struct *)ckalloc(sizeof(struct rst_struct));
115 rst->current_ff = NULL;
117 /* store object name */
118 strcpym(&rst->name, argv[1]);
120 /* create the Tcl interface to the command */
121 if(debug)printf("create object %s\n",argv[1]);
122 Tcl_CreateCommand(interp, rst->name, Chemsh_RestraintObjCmd, (ClientData) rst, (Tcl_CmdDeleteProc *) NULL);
124 /* Access initial structure */
125 GetArg(s,"coords","restraints");
126 rst->lc = get_objlist(s,"fragment",CHEMSH_OBJ_EXISTS,0);
127 if(!rst->lc){
128 printf("bad tag coords=%s\n",s);
129 return TCL_ERROR;
130 }
131 rst->c = (Frag) rst->lc->object->data;
132 rst->natms = FRAG_get_number_of_atoms(rst->c);
134 /* Gradient matrix to accumulate into (3n*1) */
135 GetArg(s,"gradient","restraints");
136 rst->lg = get_objlist(s,"matrix",CHEMSH_OBJ_EXISTS,0);
137 if(!rst->lg){
138 printf("%s: bad matrix tag %s\n", cmd,s);
139 return TCL_ERROR;
140 }
141 rst->g = (Matrix) rst->lg->object->data;
143 /* energy matrix (1*1) */
144 GetArg(s,"energy","restraints");
145 rst->le = get_objlist(s,"matrix",CHEMSH_OBJ_EXISTS,0);
146 if(!rst->le){
147 printf("%s: bad matrix tag %s\n", cmd,s);
148 return TCL_ERROR;
149 }
150 rst->e = (Matrix) rst->le->object->data;
152 /*
153 Now parse requests for terms
154 Atoms are identified by integer labels (counting from 1)
155 */
157 sprintf(tclbuff,"rst_parm(%s,%s)",rst->name,"restraints");
158 if((s = Tcl_GetVar(interp,tclbuff,0)) == NULL){
159 printf("no restraints setting in restraints module\n");
160 return TCL_ERROR;
161 }
162 if(debug)printf("restraints: %s\n",s);
164 Tcl_SplitList(interp,s,&rst->nrest,&work);
165 if(debug)printf("parsing %d restraints\n",rst->nrest);
167 for(j=0;j<rst->nrest;j++){
169 Tcl_SplitList(interp,work[j],&nfld,&work2);
171 ff=new_rst_ff(rst);
172 ff->index = j+1;
174 if(!strcmp(work2[0],"bond")){
176 /* harmonic distance restraint - expected format is bond i j r0 K */
178 if ( nfld != 5 ) {
179 printf("restraint: invalid format for restraint bond\n");
180 printf("number of field exprected: 5 (bond atom1 atom2 r0 K) \n");
181 printf("number of fields found: %d\n",nfld);
182 return TCL_ERROR;
183 }
185 if(Tcl_GetInt(interp,work2[1],&itmp)){
186 printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]);
187 return TCL_ERROR;
188 }
189 ff->i1 = itmp-1;
191 if(Tcl_GetInt(interp,work2[2],&itmp)){
192 printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]);
193 return TCL_ERROR;
194 }
195 ff->i2 = itmp-1;
197 ff->type = BOND;
199 if(Tcl_GetDouble(interp,work2[3],&atmp)){
200 printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]);
201 return TCL_ERROR;
202 }
203 ff->r0 = atmp;
205 if(Tcl_GetDouble(interp,work2[4],&atmp)){
206 printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]);
207 return TCL_ERROR;
208 }
209 ff->fc = atmp;
211 } else if(!strcmp(work2[0],"elec")){
213 /* electrostatic term - expected format is elec i j chargeprod
214 not a restraint, but added in as it may be useful for other purposes.....
215 */
217 if ( nfld != 4 ) {
218 printf("restraint: invalid format for restraint elec\n");
219 printf("number of field exprected: 4 (elec i j chargeprod) \n");
220 printf("number of fields found: %d\n",nfld);
221 return TCL_ERROR;
222 }
224 if(Tcl_GetInt(interp,work2[1],&itmp)){
225 printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]);
226 return TCL_ERROR;
227 }
228 ff->i1 = itmp-1;
230 if(Tcl_GetInt(interp,work2[2],&itmp)){
231 printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]);
232 return TCL_ERROR;
233 }
234 ff->i2 = itmp-1;
236 ff->type = ELECTROSTATIC;
238 if(Tcl_GetDouble(interp,work2[3],&atmp)){
239 printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]);
240 return TCL_ERROR;
241 }
242 ff->chargeprod = atmp;
244 }else if(!strcmp(work2[0],"angle")){
246 /* harmonic angle restraint - expected format is angle i j k r0 K */
248 if ( nfld != 6 ) {
249 printf("restraint: invalid format for restraint angle\n");
250 printf("number of field exprected: 6 (angle i j k r0 K) \n");
251 printf("number of fields found: %d\n",nfld);
252 return TCL_ERROR;
253 }
256 if(Tcl_GetInt(interp,work2[1],&itmp)){
257 printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]);
258 return TCL_ERROR;
259 }
260 ff->i1 = itmp-1;
262 if(Tcl_GetInt(interp,work2[2],&itmp)){
263 printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]);
264 return TCL_ERROR;
265 }
266 ff->i2 = itmp-1;
268 if(Tcl_GetInt(interp,work2[3],&itmp)){
269 printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]);
270 return TCL_ERROR;
271 }
272 ff->i3 = itmp-1;
273 ff->type = ANG;
274 if(Tcl_GetDouble(interp,work2[4],&atmp)){
275 printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]);
276 return TCL_ERROR;
277 }
278 /* Store internally r0 in radians */
279 #define cnv 57.29577951
280 ff->r0 = atmp / cnv;
282 if(Tcl_GetDouble(interp,work2[5],&atmp)){
283 printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]);
284 return TCL_ERROR;
285 }
286 ff->fc = atmp;
288 }else if(!strcmp(work2[0],"torsion")){
290 /* harmonic torsion restraint - expected format is torsion i j k l r0 K */
292 if ( nfld != 7 ) {
293 printf("restraint: invalid format for restraint torsion\n");
294 printf("number of field exprected: 7 (torsion i j k l r0 K) \n");
295 printf("number of fields found: %d\n",nfld);
296 return TCL_ERROR;
297 }
299 if(Tcl_GetInt(interp,work2[1],&itmp)){
300 printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]);
301 return TCL_ERROR;
302 }
303 ff->i1 = itmp-1;
305 if(Tcl_GetInt(interp,work2[2],&itmp)){
306 printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]);
307 return TCL_ERROR;
308 }
309 ff->i2 = itmp-1;
311 if(Tcl_GetInt(interp,work2[3],&itmp)){
312 printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]);
313 return TCL_ERROR;
314 }
315 ff->i3 = itmp-1;
316 if(Tcl_GetInt(interp,work2[4],&itmp)){
317 printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]);
318 return TCL_ERROR;
319 }
320 ff->i4 = itmp-1;
322 ff->type = TOR;
323 if(Tcl_GetDouble(interp,work2[5],&atmp)){
324 printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]);
325 return TCL_ERROR;
326 }
327 /* Store internally r0 in radians */
328 ff->r0 = atmp / cnv;
330 if(Tcl_GetDouble(interp,work2[6],&atmp)){
331 printf("restraint: invalid value for restraint parameter field %d (%s)\n",6,work2[6]);
332 return TCL_ERROR;
333 }
334 ff->fc = atmp;
335 }else if(!strcmp(work2[0],"bonddiff2")){
337 /* difference between two bond lengths restraint - expected format is
338 bonddiff2 i j k l r0 K */
340 if ( nfld != 7 ) {
341 printf("restraint: invalid format for restraint bonddiff2\n");
342 printf("number of field exprected: 7 (bonddiff2 i j k l r0 K) \n");
343 printf("number of fields found: %d\n",nfld);
344 return TCL_ERROR;
345 }
347 if(Tcl_GetInt(interp,work2[1],&itmp)){
348 printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]);
349 return TCL_ERROR;
350 }
351 ff->i1 = itmp-1;
352 if(Tcl_GetInt(interp,work2[2],&itmp)){
353 printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]);
354 return TCL_ERROR;
355 }
356 ff->i2 = itmp-1;
357 if(Tcl_GetInt(interp,work2[3],&itmp)){
358 printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]);
359 return TCL_ERROR;
360 }
361 ff->i3 = itmp-1;
362 if(Tcl_GetInt(interp,work2[4],&itmp)){
363 printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]);
364 return TCL_ERROR;
365 }
366 ff->i4 = itmp-1;
368 ff->type = BONDDIFF2;
370 if(Tcl_GetDouble(interp,work2[5],&atmp)){
371 printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]);
372 return TCL_ERROR;
373 }
374 ff->r0 = atmp;
376 if(Tcl_GetDouble(interp,work2[6],&atmp)){
377 printf("restraint: invalid value for restraint parameter field %d (%s)\n",6,work2[6]);
378 return TCL_ERROR;
379 }
380 ff->fc = atmp;
382 }else if(!strcmp(work2[0],"bonddiff3")){
384 /* difference between three bond lengths restraint - expected format is
385 bonddiff3 i j k l m n r0 K */
387 if ( nfld != 9 ) {
388 printf("restraint: invalid format for restraint bonddiff3\n");
389 printf("number of field exprected: 9 (bonddiff3 i j k l m n r0 K) \n");
390 printf("number of fields found: %d\n",nfld);
391 return TCL_ERROR;
392 }
394 if(Tcl_GetInt(interp,work2[1],&itmp)){
395 printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]);
396 return TCL_ERROR;
397 }
398 ff->i1 = itmp-1;
399 if(Tcl_GetInt(interp,work2[2],&itmp)){
400 printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]);
401 return TCL_ERROR;
402 }
403 ff->i2 = itmp-1;
404 if(Tcl_GetInt(interp,work2[3],&itmp)){
405 printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]);
406 return TCL_ERROR;
407 }
408 ff->i3 = itmp-1;
409 if(Tcl_GetInt(interp,work2[4],&itmp)){
410 printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]);
411 return TCL_ERROR;
412 }
413 ff->i4 = itmp-1;
414 if(Tcl_GetInt(interp,work2[5],&itmp)){
415 printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]);
416 return TCL_ERROR;
417 }
418 ff->i5 = itmp-1;
419 if(Tcl_GetInt(interp,work2[6],&itmp)){
420 printf("restraint: invalid value for restraint parameter field %d (%s)\n",6,work2[6]);
421 return TCL_ERROR;
422 }
423 ff->i6 = itmp-1;
425 ff->type = BONDDIFF3;
427 if(Tcl_GetDouble(interp,work2[7],&atmp)){
428 printf("restraint: invalid value for restraint parameter field %d (%s)\n",7,work2[7]);
429 return TCL_ERROR;
430 }
431 ff->r0 = atmp;
433 if(Tcl_GetDouble(interp,work2[8],&atmp)){
434 printf("restraint: invalid value for restraint parameter field %d (%s)\n",8,work2[8]);
435 return TCL_ERROR;
436 }
437 ff->fc = atmp;
439 }else if(!strcmp(work2[0],"bonddiff4")){
441 /* difference between four bond lengths restraint - expected format is
442 bonddiff3 i j k l m n o p r0 K
443 E=K/2 * ( d(ij)+d(kl)-d(mn)-d(op) - r0 )**2 */
445 if ( nfld != 11 ) {
446 printf("restraint: invalid format for restraint bonddiff4\n");
447 printf("number of field exprected: 11 (bonddiff3 i j k l m n o p r0 K) \n");
448 printf("number of fields found: %d\n",nfld);
449 return TCL_ERROR;
450 }
452 if(Tcl_GetInt(interp,work2[1],&itmp)){
453 printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]);
454 return TCL_ERROR;
455 }
456 ff->i1 = itmp-1;
457 if(Tcl_GetInt(interp,work2[2],&itmp)){
458 printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]);
459 return TCL_ERROR;
460 }
461 ff->i2 = itmp-1;
462 if(Tcl_GetInt(interp,work2[3],&itmp)){
463 printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]);
464 return TCL_ERROR;
465 }
466 ff->i3 = itmp-1;
467 if(Tcl_GetInt(interp,work2[4],&itmp)){
468 printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]);
469 return TCL_ERROR;
470 }
471 ff->i4 = itmp-1;
472 if(Tcl_GetInt(interp,work2[5],&itmp)){
473 printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]);
474 return TCL_ERROR;
475 }
476 ff->i5 = itmp-1;
477 if(Tcl_GetInt(interp,work2[6],&itmp)){
478 printf("restraint: invalid value for restraint parameter field %d (%s)\n",6,work2[6]);
479 return TCL_ERROR;
480 }
481 ff->i6 = itmp-1;
482 if(Tcl_GetInt(interp,work2[7],&itmp)){
483 printf("restraint: invalid value for restraint parameter field %d (%s)\n",7,work2[7]);
484 return TCL_ERROR;
485 }
486 ff->i7 = itmp-1;
487 if(Tcl_GetInt(interp,work2[8],&itmp)){
488 printf("restraint: invalid value for restraint parameter field %d (%s)\n",8,work2[8]);
489 return TCL_ERROR;
490 }
491 ff->i8 = itmp-1;
493 ff->type = BONDDIFF4;
495 if(Tcl_GetDouble(interp,work2[9],&atmp)){
496 printf("restraint: invalid value for restraint parameter field %d (%s)\n",7,work2[9]);
497 return TCL_ERROR;
498 }
499 ff->r0 = atmp;
501 if(Tcl_GetDouble(interp,work2[10],&atmp)){
502 printf("restraint: invalid value for restraint parameter field %d (%s)\n",8,work2[10]);
503 return TCL_ERROR;
504 }
505 ff->fc = atmp;
507 }else if(!strcmp(work2[0],"sbc_atom")){
509 /* Spherical boundary conditions around an atom - expected format is
510 sbc_atom i r0 K E= Sum ( 0.5 * K * (| r_a - r_i | - r0 )^2 ) */
512 if ( nfld != 4 ) {
513 printf("restraint: invalid format for restraint sbc_atom\n");
514 printf("number of field exprected: 4 (sbc_atom i r0 K) \n");
515 printf("number of fields found: %d\n",nfld);
516 return TCL_ERROR;
517 }
519 if(Tcl_GetInt(interp,work2[1],&itmp)){
520 printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]);
521 return TCL_ERROR;
522 }
523 ff->i1 = itmp-1;
524 if(Tcl_GetDouble(interp,work2[2],&atmp)){
525 printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]);
526 return TCL_ERROR;
527 }
528 ff->r0 = atmp;
530 if(Tcl_GetDouble(interp,work2[3],&atmp)){
531 printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]);
532 return TCL_ERROR;
533 }
534 ff->fc = atmp;
535 ff->type = SBC_ATOM;
537 }else if(!strcmp(work2[0],"sbc")){
539 /* Spherical boundary conditions around an geometric center c,
540 specified by x,y,z - expected format is sbc x y z r0 K
541 E= Sum ( 0.5 * K * (| r_a - c | - r0 )^2 ) */
543 if ( nfld != 6 ) {
544 printf("restraint: invalid format for restraint sbc_atom\n");
545 printf("number of field exprected: 6 (sbc x y z r0 K) \n");
546 printf("number of fields found: %d\n",nfld);
547 return TCL_ERROR;
548 }
550 if(Tcl_GetDouble(interp,work2[1],&atmp)){
551 printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]);
552 return TCL_ERROR;
553 }
554 ff->x = atmp;
555 if(Tcl_GetDouble(interp,work2[2],&atmp)){
556 printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]);
557 return TCL_ERROR;
558 }
559 ff->y = atmp;
560 if(Tcl_GetDouble(interp,work2[3],&atmp)){
561 printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]);
562 return TCL_ERROR;
563 }
564 ff->z = atmp;
565 if(Tcl_GetDouble(interp,work2[4],&atmp)){
566 printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]);
567 return TCL_ERROR;
568 }
569 ff->r0 = atmp;
571 if(Tcl_GetDouble(interp,work2[5],&atmp)){
572 printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]);
573 return TCL_ERROR;
574 }
575 ff->fc = atmp;
576 ff->type = SBC;
578 }else if(!strcmp(work2[0],"sbc_individual")){
580 /* Spherical boundary conditions around an geometric center c for an explicit atom,
581 specified by x,y,z - expected format is sbc_individual atom_number x y z r0 K
582 E= Sum ( 0.5 * K * (| r_a - c | - r0 )^2 ) */
584 if ( nfld != 7 ) {
585 printf("restraint: invalid format for restraint sbc_individual\n");
586 printf("number of field exprected: 7 (sbc_individual atom_number x y z r0 K) \n");
587 printf("number of fields found: %d\n",nfld);
588 return TCL_ERROR;
589 }
591 if(Tcl_GetInt(interp,work2[1],&itmp)){
592 printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]);
593 return TCL_ERROR;
594 }
596 ff->i1 = itmp-1;
598 if(Tcl_GetDouble(interp,work2[2],&atmp)){
599 printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]);
600 return TCL_ERROR;
601 }
602 ff->x = atmp;
603 if(Tcl_GetDouble(interp,work2[3],&atmp)){
604 printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]);
605 return TCL_ERROR;
606 }
607 ff->y = atmp;
608 if(Tcl_GetDouble(interp,work2[4],&atmp)){
609 printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]);
610 return TCL_ERROR;
611 }
612 ff->z = atmp;
613 if(Tcl_GetDouble(interp,work2[5],&atmp)){
614 printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]);
615 return TCL_ERROR;
616 }
617 ff->r0 = atmp;
619 if(Tcl_GetDouble(interp,work2[6],&atmp)){
620 printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]);
621 return TCL_ERROR;
622 }
623 ff->fc = atmp;
624 ff->type = SBC_IND;
626 }else if(!strcmp(work2[0],"cart_atom")){
628 /* Restrains an atom around its original position
629 cart_atom iat K
630 E= 0.5 * K * (| r_a - c | )^2 */
632 if ( nfld != 3 ) {
633 printf("restraint: invalid format for restraint cart_atom\n");
634 printf("number of field exprected: 3 (cart_atom atom_number K) \n");
635 printf("number of fields found: %d\n",nfld);
636 return TCL_ERROR;
637 }
639 if(Tcl_GetInt(interp,work2[1],&itmp)){
640 printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]);
641 return TCL_ERROR;
642 }
643 ff->i1 = itmp-1;
644 ff->x = rst->c->atoms[itmp-1].pos.x[0];
645 ff->y = rst->c->atoms[itmp-1].pos.x[1];
646 ff->z = rst->c->atoms[itmp-1].pos.x[2];
648 ff->r0 = 0.;
650 if(Tcl_GetDouble(interp,work2[2],&atmp)){
651 printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]);
652 return TCL_ERROR;
653 }
654 ff->fc = atmp;
655 ff->type = CART_ATOM;
657 }else if(!strcmp(work2[0],"cart")){
659 /* Restrains an atom around the position x y z
660 cart iat x y z K
661 E= 0.5 * K * (| r_a - c | )^2 */
663 if ( nfld != 6 ) {
664 printf("restraint: invalid format for restraint cart\n");
665 printf("number of field exprected: 6 (cart atom_number x y z K) \n");
666 printf("number of fields found: %d\n",nfld);
667 return TCL_ERROR;
668 }
670 if(Tcl_GetInt(interp,work2[1],&itmp)){
671 printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]);
672 return TCL_ERROR;
673 }
674 ff->i1 = itmp-1;
676 if(Tcl_GetDouble(interp,work2[2],&atmp)){
677 printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]);
678 return TCL_ERROR;
679 }
680 ff->x = atmp;
681 if(Tcl_GetDouble(interp,work2[3],&atmp)){
682 printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]);
683 return TCL_ERROR;
684 }
685 ff->y = atmp;
686 if(Tcl_GetDouble(interp,work2[4],&atmp)){
687 printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]);
688 return TCL_ERROR;
689 }
690 ff->z = atmp;
692 ff->r0 = 0.;
694 if(Tcl_GetDouble(interp,work2[5],&atmp)){
695 printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]);
696 return TCL_ERROR;
697 }
698 ff->fc = atmp;
699 ff->type = CART;
701 } else {
702 printf("restraint: invalid value for restraint type (%s)\n",work2[0]);
703 return TCL_ERROR;
704 }
705 ckfree((char *) work2);
706 }
707 ckfree((char *) work);
709 printf("Restraints Atoms (object name %s) :\n",rst->name);
710 for(ff=rst->last_ff;ff;ff=ff->prev){
711 switch(ff->type){
712 case BOND:
713 printf("Harmonic Distance %5d %5d r0 = %f fc = %f \n",
714 ff->i1+1,
715 ff->i2+1,
716 ff->r0,
717 ff->fc );
718 break;
719 case ANG:
720 printf("Harmonic Angle %5d %5d %5d r0 = %f fc = %f \n",
721 ff->i1+1,
722 ff->i2+1,
723 ff->i3+1,
724 ff->r0,
725 ff->fc );
726 break;
727 case ELECTROSTATIC:
728 printf("Electrostatic %5d %5d qiqj = %f \n",
729 ff->i1+1,
730 ff->i2+1,
731 ff->chargeprod);
732 break;
733 case TOR:
734 printf("Torsion %5d %5d %5d %5d r0 = %f fc = %f \n",
735 ff->i1+1,
736 ff->i2+1,
737 ff->i3+1,
738 ff->i4+1,
739 ff->r0,
740 ff->fc );
741 break;
742 case BONDDIFF2:
743 printf("Bonddiff-2 %5d %5d %5d %5d r0 = %f fc = %f \n",
744 ff->i1+1,
745 ff->i2+1,
746 ff->i3+1,
747 ff->i4+1,
748 ff->r0,
749 ff->fc );
750 break;
751 case BONDDIFF3:
752 printf("Bonddiff-3 %5d %5d %5d %5d %5d %5d r0 = %f fc = %f \n",
753 ff->i1+1,
754 ff->i2+1,
755 ff->i3+1,
756 ff->i4+1,
757 ff->i5+1,
758 ff->i6+1,
759 ff->r0,
760 ff->fc );
761 break;
762 case BONDDIFF4:
763 printf("Bonddiff-4 %5d %5d %5d %5d %5d %5d %5d %5d r0 = %f fc = %f \n",
764 ff->i1+1,
765 ff->i2+1,
766 ff->i3+1,
767 ff->i4+1,
768 ff->i5+1,
769 ff->i6+1,
770 ff->i7+1,
771 ff->i8+1,
772 ff->r0,
773 ff->fc );
774 break;
775 case SBC_ATOM:
776 printf("SBC Atom %5d r0 = %f fc = %f \n",
777 ff->i1+1,
778 ff->r0,
779 ff->fc );
780 break;
781 case SBC:
782 printf("SBC %f %f %f r0 = %f fc = %f \n",
783 ff->x,
784 ff->y,
785 ff->z,
786 ff->r0,
787 ff->fc );
788 break;
789 case SBC_IND:
790 printf("SBC_individual %5d %f %f %f r0 = %f fc = %f \n",
791 ff->i1+1,
792 ff->x,
793 ff->y,
794 ff->z,
795 ff->r0,
796 ff->fc );
797 break;
798 case CART_ATOM:
799 printf("Cart Atom %5d %10.5f %10.5f %10.5f fc = %f \n",
800 ff->i1+1,
801 ff->x,
802 ff->y,
803 ff->z,
804 ff->fc );
805 break;
806 case CART:
807 printf("Cart %5d %10.5f %10.5f %10.5f fc = %f \n",
808 ff->i1+1,
809 ff->x,
810 ff->y,
811 ff->z,
812 ff->fc );
813 break;
815 }
816 }
817 return TCL_OK;
818 }
820 /*
821 Implements the commands associated with an restraint object
822 The first argument is the required method
823 */
824 int
825 Chemsh_RestraintObjCmd(ClientData clientData,
826 Tcl_Interp *interp,
827 int argc,
828 const char *argv[])
829 {
831 /* local variables */
832 char tclbuff[STRBUFFLEN];
833 char strerr[5];
834 char c;
835 int i, j, length;
836 struct rst_struct *rst;
837 Ff ff;
838 int itmp;
840 rst = (struct rst_struct *) clientData;
842 /* default Tcl result */
843 strcpy(strerr,"0");
845 if (argc == 1){
846 printf("optimiser object %s invoked without a command\n",rst->name);
847 return TCL_ERROR;
848 }
850 c = argv[1][0]; /* First letter of command */
851 length = strlen(argv[1]);
853 sprintf(tclbuff,"global rst_parm");
854 if (Tcl_Eval(interp,tclbuff) != TCL_OK) return TCL_ERROR;
856 if(debug) printf("%s Invoking method: %s\n",rst->name, argv[1]);
858 if ((c == 'd') && ( strncmp(argv[1],"destroy", length) == 0)) {
860 /* nothing to do really... maybe release rst structure ??? */
861 rel_objlist(rst->lc);
862 rel_objlist(rst->lg);
863 rel_objlist(rst->le);
865 /* delete_restraint_memory(rst); */
866 /*pop_module_name(interp); */
868 } else if ((c == 'c') && ( strncmp(argv[1],"configure", length) == 0)) {
870 /* Configure */
871 if(rst_newconfig(interp,rst,argc-2,argv+2)) {
872 return TCL_ERROR;
873 }
874 /*rst_print_parms(rst);*/
876 } else if ((c == 'o') && ( strncmp(argv[1],"output", length) == 0)) {
878 /* Print current values of restraint values and energies and force terms */
880 } else if ((c == 'a') && ( strncmp(argv[1],"accumulate", length) == 0)) {
882 /* Accumulate and print results */
883 for(ff=rst->last_ff;ff;ff=ff->prev){
884 calc_rst_ff(ff, rst->c, rst->e, rst->g);
885 }
887 } else {
889 printf("unrecognised method for restraint object %s : %s\n",rst->name, argv[1]);
890 sprintf(tclbuff,"unrecognised method for restraint object %s : %s\n",rst->name, argv[1]);
891 Tcl_SetResult(interp, tclbuff, TCL_VOLATILE);
892 return TCL_ERROR;
894 }
896 /* clean exits */
897 Tcl_SetResult(interp, strerr, TCL_VOLATILE);
898 return TCL_OK;
900 }
902 /* Process configuration options */
904 int rst_newconfig(Tcl_Interp *interp,
905 struct rst_struct *rst,
906 int argc,
907 const char **argv) {
909 char tclbuff[STRBUFFLEN];
910 int i;
912 for(i=0;i<argc;i++){
913 Tcl_SetVar(interp,"argtmp",argv[i],TCL_APPEND_VALUE | TCL_LIST_ELEMENT);
914 }
916 sprintf(tclbuff,"rst_parse %s {temperature pressure timestep ensemble trajectory_type taut taup energy_unit} $argtmp ",rst->name);
918 if(Tcl_Eval(interp,tclbuff)){
919 printf("%s: configure failed - %s\n",rst->name,interp->result);
920 return TCL_ERROR;
921 }
922 rst_set_config(interp,rst);
924 printf("%s: configure completed: ",rst->name);
925 for(i=0;i<argc;i++){
926 printf(" ");
927 printf(argv[i]);
928 }
929 printf("\n");
931 return TCL_OK;
932 }
935 /* Set configurable parameters */
937 int
938 rst_set_config(Tcl_Interp *interp,
939 struct rst_struct *rst) {
941 /* nothing here yet */
943 return 0;
944 }
946 int
947 calc_rst_ff(Ff ff,
948 Frag f,
949 Matrix e,
950 Matrix g)
951 {
952 double d, ee, r12, r32, r34, r56, r78, phi, cosphi, sinphi;
953 double torpos[4][3],torval,torder[4][3],fact,pi,pressure;
954 int icor,iat;
955 struct vector_struct v12, v32, f12, f1, f3;
956 struct atom_struct *a1, *a2, *a3, *a4, *a5, *a6, *a7, *a8,b1,b2;
957 struct vector_struct v34, f34, v56, v78, f56, f78;
959 switch(ff->type){
960 case BOND:
961 a1 = &f->atoms[ff->i1];
962 a2 = &f->atoms[ff->i2];
963 r12 = vdist(a1->pos,a2->pos);
964 v12 = vnorm(vdiff(a2->pos,a1->pos));
965 d = r12 - ff->r0;
966 f12 = vscal(v12,-d*ff->fc);
967 if(g->d){
968 g->d[ff->i1*3 + 0] += f12.x[0];
969 g->d[ff->i1*3 + 1] += f12.x[1];
970 g->d[ff->i1*3 + 2] += f12.x[2];
971 g->d[ff->i2*3 + 0] -= f12.x[0];
972 g->d[ff->i2*3 + 1] -= f12.x[1];
973 g->d[ff->i2*3 + 2] -= f12.x[2];
974 }
975 ee = 0.5*d*d*ff->fc;
976 e->d[0] += ee;
977 if(debug)printf("r,energy in calc_ff %lf %10.4e\n",r12,ee);
978 printf("Restraint %2d dist : value %11.4e energy %11.4e force %11.4e\n",ff->index,r12,ee,-d*ff->fc);
979 break;
981 case ELECTROSTATIC:
982 a1 = &f->atoms[ff->i1];
983 a2 = &f->atoms[ff->i2];
984 r12 = vdist(a1->pos,a2->pos);
985 v12 = vnorm(vdiff(a2->pos,a1->pos));
986 f12 = vscal(v12,ff->chargeprod*(1.0/(r12*r12)));
987 if(g->d){
988 g->d[ff->i1*3 + 0] += f12.x[0];
989 g->d[ff->i1*3 + 1] += f12.x[1];
990 g->d[ff->i1*3 + 2] += f12.x[2];
991 g->d[ff->i2*3 + 0] -= f12.x[0];
992 g->d[ff->i2*3 + 1] -= f12.x[1];
993 g->d[ff->i2*3 + 2] -= f12.x[2];
994 }
995 ee = ff->chargeprod*(1.0/r12);
996 e->d[0] += ee;
997 if(debug)printf("chargeprod, r,energy in calc_ff %f %lf %10.4e\n",ff->chargeprod,r12,ee);
998 printf("Restraint %2d elec : value %11.4e energy %11.4e\n",ff->index,r12,ee);
999 break;
1001 case ANG:
1002 a1 = &f->atoms[ff->i1];
1003 a2 = &f->atoms[ff->i2];
1004 a3 = &f->atoms[ff->i3];
1005 r12 = vdist(a1->pos,a2->pos);
1006 r32 = vdist(a3->pos,a2->pos);
1007 v12 = vnorm(vdiff(a2->pos,a1->pos));
1008 v32 = vnorm(vdiff(a2->pos,a3->pos));
1009 cosphi = vdot(v12,v32);
1010 phi = acos(cosphi);
1011 sinphi = sin(phi);
1012 d = phi - ff->r0;
1013 f1 = vscal(vdiff(vscal(v12,cosphi),v32), d*ff->fc* 1.0/(r12*sinphi));
1014 f3 = vscal(vdiff(vscal(v32,cosphi),v12), d*ff->fc* 1.0/(r32*sinphi));
1015 if(debug){
1016 vprint("v12",v12);
1017 vprint("v32",v32);
1018 vprint("f1",f1);
1019 vprint("f3",f3);
1020 }
1021 if(g->d){
1022 g->d[ff->i1*3 + 0] -= f1.x[0];
1023 g->d[ff->i1*3 + 1] -= f1.x[1];
1024 g->d[ff->i1*3 + 2] -= f1.x[2];
1025 g->d[ff->i3*3 + 0] -= f3.x[0];
1026 g->d[ff->i3*3 + 1] -= f3.x[1];
1027 g->d[ff->i3*3 + 2] -= f3.x[2];
1028 g->d[ff->i2*3 + 0] += (f1.x[0] + f3.x[0]);
1029 g->d[ff->i2*3 + 1] += (f1.x[1] + f3.x[1]);
1030 g->d[ff->i2*3 + 2] += (f1.x[2] + f3.x[2]);
1031 }
1032 ee = 0.5 * d * d * ff->fc;
1033 if(debug)printf("phi,energy in calc_ff %lf %10.4e\n",phi,ee);
1034 printf("Restraint %2d angl : value %11.4e energy %11.4e force %11.4e\n",ff->index,phi,ee,-d*ff->fc);
1035 e->d[0] += ee;
1036 break;
1037 case TOR:
1038 /* map to fortran routines in the dynamic-module */
1039 for(icor=1;icor<=3;icor++) {
1040 a1 = &f->atoms[ff->i1];
1041 torpos[0][icor-1]=a1->pos.x[icor-1];
1042 a1 = &f->atoms[ff->i2];
1043 torpos[1][icor-1]=a1->pos.x[icor-1];
1044 a1 = &f->atoms[ff->i3];
1045 torpos[2][icor-1]=a1->pos.x[icor-1];
1046 a1 = &f->atoms[ff->i4];
1047 torpos[3][icor-1]=a1->pos.x[icor-1];
1048 }
1049 if(g->d){
1050 torsionderivative_(torpos,&phi,torder);
1051 } else {
1052 torsionval_(torpos,&phi);
1053 }
1054 d = phi - ff->r0;
1055 /* check for a 360 deg rotation and avoid it */
1056 pi = 4. * atan(1.);
1057 if( d > 7. || d < -7.) {
1058 printf("Error: torsion values differ too much!\n");
1059 printf("Current value %f, centre of the restraint %f\n",phi,ff->r0);
1060 exit(1);
1061 }
1062 if(d > 3.1415 ) d=d-2.*pi;
1063 if(d < -3.1415 ) d=d+2.*pi;
1065 if(g->d){
1066 fact= ff->fc * d;
1067 g->d[ff->i1*3 + 0] += fact*torder[0][0];
1068 g->d[ff->i1*3 + 1] += fact*torder[0][1];
1069 g->d[ff->i1*3 + 2] += fact*torder[0][2];
1070 g->d[ff->i2*3 + 0] += fact*torder[1][0];
1071 g->d[ff->i2*3 + 1] += fact*torder[1][1];
1072 g->d[ff->i2*3 + 2] += fact*torder[1][2];
1073 g->d[ff->i3*3 + 0] += fact*torder[2][0];
1074 g->d[ff->i3*3 + 1] += fact*torder[2][1];
1075 g->d[ff->i3*3 + 2] += fact*torder[2][2];
1076 g->d[ff->i4*3 + 0] += fact*torder[3][0];
1077 g->d[ff->i4*3 + 1] += fact*torder[3][1];
1078 g->d[ff->i4*3 + 2] += fact*torder[3][2];
1079 }
1080 ee = 0.5 * d * d * ff->fc;
1081 if(debug)printf("Torsion: phi, energy in calc_ff %lf %10.4e\n",phi,ee);
1082 printf("Restraint %2d tors : value %11.4e energy %11.4e force %11.4e\n",ff->index,phi,ee,-d*ff->fc);
1083 e->d[0] += ee;
1084 break;
1085 case BONDDIFF2:
1086 a1 = &f->atoms[ff->i1];
1087 a2 = &f->atoms[ff->i2];
1088 a3 = &f->atoms[ff->i3];
1089 a4 = &f->atoms[ff->i4];
1090 r12 = vdist(a1->pos,a2->pos);
1091 r34 = vdist(a3->pos,a4->pos);
1092 v12 = vnorm(vdiff(a2->pos,a1->pos));
1093 v34 = vnorm(vdiff(a4->pos,a3->pos));
1094 d = r12 - r34 - ff->r0;
1095 f12 = vscal(v12,-d*ff->fc);
1096 f34 = vscal(v34,-d*ff->fc);
1097 if(g->d){
1098 g->d[ff->i1*3 + 0] += f12.x[0];
1099 g->d[ff->i1*3 + 1] += f12.x[1];
1100 g->d[ff->i1*3 + 2] += f12.x[2];
1101 g->d[ff->i2*3 + 0] -= f12.x[0];
1102 g->d[ff->i2*3 + 1] -= f12.x[1];
1103 g->d[ff->i2*3 + 2] -= f12.x[2];
1104 g->d[ff->i3*3 + 0] -= f34.x[0];
1105 g->d[ff->i3*3 + 1] -= f34.x[1];
1106 g->d[ff->i3*3 + 2] -= f34.x[2];
1107 g->d[ff->i4*3 + 0] += f34.x[0];
1108 g->d[ff->i4*3 + 1] += f34.x[1];
1109 g->d[ff->i4*3 + 2] += f34.x[2];
1110 }
1111 ee = 0.5*d*d*ff->fc;
1112 e->d[0] += ee;
1113 printf("Restraint %2d bd2 : value %11.4e energy %11.4e force %11.4e\n",ff->index,r12-r34,ee,-d*ff->fc);
1114 break;
1116 case BONDDIFF3:
1117 a1 = &f->atoms[ff->i1];
1118 a2 = &f->atoms[ff->i2];
1119 a3 = &f->atoms[ff->i3];
1120 a4 = &f->atoms[ff->i4];
1121 a5 = &f->atoms[ff->i5];
1122 a6 = &f->atoms[ff->i6];
1123 r12 = vdist(a1->pos,a2->pos);
1124 r34 = vdist(a3->pos,a4->pos);
1125 r56 = vdist(a5->pos,a6->pos);
1126 v12 = vnorm(vdiff(a2->pos,a1->pos));
1127 v34 = vnorm(vdiff(a4->pos,a3->pos));
1128 v56 = vnorm(vdiff(a6->pos,a5->pos));
1129 d = r12 - r34 -r56 - ff->r0;
1130 f12 = vscal(v12,-d*ff->fc);
1131 f34 = vscal(v34,-d*ff->fc);
1132 f56 = vscal(v56,-d*ff->fc);
1133 if(g->d){
1134 g->d[ff->i1*3 + 0] += f12.x[0];
1135 g->d[ff->i1*3 + 1] += f12.x[1];
1136 g->d[ff->i1*3 + 2] += f12.x[2];
1137 g->d[ff->i2*3 + 0] -= f12.x[0];
1138 g->d[ff->i2*3 + 1] -= f12.x[1];
1139 g->d[ff->i2*3 + 2] -= f12.x[2];
1140 g->d[ff->i3*3 + 0] -= f34.x[0];
1141 g->d[ff->i3*3 + 1] -= f34.x[1];
1142 g->d[ff->i3*3 + 2] -= f34.x[2];
1143 g->d[ff->i4*3 + 0] += f34.x[0];
1144 g->d[ff->i4*3 + 1] += f34.x[1];
1145 g->d[ff->i4*3 + 2] += f34.x[2];
1146 g->d[ff->i5*3 + 0] -= f56.x[0];
1147 g->d[ff->i5*3 + 1] -= f56.x[1];
1148 g->d[ff->i5*3 + 2] -= f56.x[2];
1149 g->d[ff->i6*3 + 0] += f56.x[0];
1150 g->d[ff->i6*3 + 1] += f56.x[1];
1151 g->d[ff->i6*3 + 2] += f56.x[2];
1152 }
1153 ee = 0.5*d*d*ff->fc;
1154 e->d[0] += ee;
1155 printf("Restraint %2d bd3 : value %11.4e energy %11.4e force %11.4e\n",ff->index,r12-r34-r56,ee,-d*ff->fc);
1156 break;
1157 case BONDDIFF4:
1158 a1 = &f->atoms[ff->i1];
1159 a2 = &f->atoms[ff->i2];
1160 a3 = &f->atoms[ff->i3];
1161 a4 = &f->atoms[ff->i4];
1162 a5 = &f->atoms[ff->i5];
1163 a6 = &f->atoms[ff->i6];
1164 a7 = &f->atoms[ff->i7];
1165 a8 = &f->atoms[ff->i8];
1166 r12 = vdist(a1->pos,a2->pos);
1167 r34 = vdist(a3->pos,a4->pos);
1168 r56 = vdist(a5->pos,a6->pos);
1169 r78 = vdist(a7->pos,a8->pos);
1170 v12 = vnorm(vdiff(a2->pos,a1->pos));
1171 v34 = vnorm(vdiff(a4->pos,a3->pos));
1172 v56 = vnorm(vdiff(a6->pos,a5->pos));
1173 v78 = vnorm(vdiff(a8->pos,a7->pos));
1174 d = r12 + r34 - r56 - r78 - ff->r0;
1175 f12 = vscal(v12,-d*ff->fc);
1176 f34 = vscal(v34,-d*ff->fc);
1177 f56 = vscal(v56,-d*ff->fc);
1178 f78 = vscal(v78,-d*ff->fc);
1179 if(g->d){
1180 g->d[ff->i1*3 + 0] += f12.x[0];
1181 g->d[ff->i1*3 + 1] += f12.x[1];
1182 g->d[ff->i1*3 + 2] += f12.x[2];
1183 g->d[ff->i2*3 + 0] -= f12.x[0];
1184 g->d[ff->i2*3 + 1] -= f12.x[1];
1185 g->d[ff->i2*3 + 2] -= f12.x[2];
1187 g->d[ff->i3*3 + 0] += f34.x[0];
1188 g->d[ff->i3*3 + 1] += f34.x[1];
1189 g->d[ff->i3*3 + 2] += f34.x[2];
1190 g->d[ff->i4*3 + 0] -= f34.x[0];
1191 g->d[ff->i4*3 + 1] -= f34.x[1];
1192 g->d[ff->i4*3 + 2] -= f34.x[2];
1194 g->d[ff->i5*3 + 0] -= f56.x[0];
1195 g->d[ff->i5*3 + 1] -= f56.x[1];
1196 g->d[ff->i5*3 + 2] -= f56.x[2];
1197 g->d[ff->i6*3 + 0] += f56.x[0];
1198 g->d[ff->i6*3 + 1] += f56.x[1];
1199 g->d[ff->i6*3 + 2] += f56.x[2];
1201 g->d[ff->i7*3 + 0] -= f78.x[0];
1202 g->d[ff->i7*3 + 1] -= f78.x[1];
1203 g->d[ff->i7*3 + 2] -= f78.x[2];
1204 g->d[ff->i8*3 + 0] += f78.x[0];
1205 g->d[ff->i8*3 + 1] += f78.x[1];
1206 g->d[ff->i8*3 + 2] += f78.x[2];
1207 }
1208 ee = 0.5*d*d*ff->fc;
1209 e->d[0] += ee;
1210 printf("Restraint %2d bd4 : value %11.4e energy %11.4e force %11.4e\n",ff->index,r12+r34-r56-r78,ee,-d*ff->fc);
1211 break;
1212 case SBC_ATOM:
1213 ee= 0.;
1214 r32=0.;
1215 for(iat=0; iat < f->natom; iat++){
1216 a1 = &f->atoms[ff->i1];
1217 a2 = &f->atoms[iat]; /* Centre */
1218 r12 = vdist(a1->pos,a2->pos);
1219 if(r12 > r32) r32=r12;
1220 d = r12 - ff->r0;
1221 if(d > 0.0) {
1222 v12 = vnorm(vdiff(a2->pos,a1->pos));
1223 f12 = vscal(v12,d*ff->fc);
1224 if(g->d){
1225 g->d[iat*3 + 0] += f12.x[0];
1226 g->d[iat*3 + 1] += f12.x[1];
1227 g->d[iat*3 + 2] += f12.x[2];
1228 g->d[ff->i1 + 0] -= f12.x[0];
1229 g->d[ff->i1 + 1] -= f12.x[1];
1230 g->d[ff->i1 + 2] -= f12.x[2];
1231 }
1232 ee += 0.5*d*d*ff->fc;
1233 }
1234 }
1235 e->d[0] += ee;
1236 printf("Restraint %2d SBC-A: value %11.4e energy %11.4e\n",ff->index,r32,ee);
1237 break;
1239 case SBC:
1240 ee= 0.;
1241 r32=0.;
1242 pressure=0.;
1243 pi = 4. * atan(1.);
1244 for(iat=0; iat < f->natom; iat++){
1245 b1 = f->atoms[iat];
1246 b2.pos.x[0]= ff->x;
1247 b2.pos.x[1]= ff->y;
1248 b2.pos.x[2]= ff->z;
1249 r12 = vdist(b1.pos,b2.pos);
1250 if(r12 > r32) r32=r12;
1251 d = r12 - ff->r0;
1252 /* ignore the center atom */
1253 if(d > 0.0) {
1254 v12 = vnorm(vdiff(b2.pos,b1.pos));
1255 f12 = vscal(v12,d*ff->fc);
1256 if(g->d){
1257 g->d[iat*3 + 0] -= f12.x[0];
1258 g->d[iat*3 + 1] -= f12.x[1];
1259 g->d[iat*3 + 2] -= f12.x[2];
1260 }
1261 pressure += d * ff->fc / (4.*pi*r12*r12);
1262 ee += 0.5*d*d*ff->fc;
1263 }
1264 }
1265 e->d[0] += ee;
1266 printf("Restraint %2d SBC : value %11.4e energy %11.4e\n",ff->index,r32,ee);
1267 /*
1268 printf("SBC pressure: %11.4e a.u. = %11.4e bar\n",pressure, pressure * 2.9421011987E+8 );
1269 */
1270 break;
1272 case SBC_IND:
1273 ee= 0.;
1274 r32=0.;
1275 pressure=0.;
1276 pi = 4. * atan(1.);
1277 b1 = f->atoms[ff->i1];
1278 b2.pos.x[0]= ff->x;
1279 b2.pos.x[1]= ff->y;
1280 b2.pos.x[2]= ff->z;
1281 r12 = vdist(b1.pos,b2.pos);
1282 if(r12 > r32) r32=r12;
1283 d = r12 - ff->r0;
1284 /* ignore the center atom */
1285 if(d > 0.0) {
1286 v12 = vnorm(vdiff(b2.pos,b1.pos));
1287 f12 = vscal(v12,d*ff->fc);
1288 if(g->d){
1289 g->d[ff->i1*3 + 0] -= f12.x[0];
1290 g->d[ff->i1*3 + 1] -= f12.x[1];
1291 g->d[ff->i1*3 + 2] -= f12.x[2];
1292 }
1293 pressure += d * ff->fc / (4.*pi*r12*r12);
1294 ee = 0.5*d*d*ff->fc;
1295 }
1296 e->d[0] += ee;
1297 if(debug)printf("Restraint %2d SBC_individual : value %11.4e energy %11.4e\n",ff->index,r32,ee);
1298 break;
1300 case CART:
1301 /* handle cart and cart_atom on an equal basis */
1302 case CART_ATOM:
1303 b1 = f->atoms[ff->i1];
1304 b2.pos.x[0]= ff->x;
1305 b2.pos.x[1]= ff->y;
1306 b2.pos.x[2]= ff->z;
1307 r12 = vdist(b1.pos,b2.pos);
1308 d = r12;
1309 if(r12 > 0. ) {
1310 v12 = vnorm(vdiff(b2.pos,b1.pos));
1311 f12 = vscal(v12,-d*ff->fc);
1312 if(g->d){
1313 g->d[ff->i1*3 + 0] += f12.x[0];
1314 g->d[ff->i1*3 + 1] += f12.x[1];
1315 g->d[ff->i1*3 + 2] += f12.x[2];
1316 }
1317 }
1318 ee = 0.5*d*d*ff->fc;
1319 e->d[0] += ee;
1320 if(debug)printf("r,energy in calc_ff %lf %10.4e\n",r12,ee);
1321 printf("Restraint %2d cart : value %11.4e energy %11.4e force %11.4e\n",ff->index,r12,ee,-d*ff->fc);
1323 break;
1325 default:
1326 printf("energy term not implemented\n");
1327 exit(1);
1328 }
1329 fflush(stdout);
1330 return 0;
1331 }
1333 /* new_ff structure, to add an ff maintaining order */
1334 Ff new_rst_ff(struct rst_struct *rst)
1335 {
1336 Ff ff;
1337 ff = New(ff_s);
1338 if(rst->current_ff){
1339 rst->current_ff->prev = ff;
1340 }else{
1341 rst->last_ff = ff;
1342 }
1343 rst->current_ff = ff;
1344 ff->prev = NULL;
1345 return ff;
1346 }
1348 /* Here starts the Block for the dispersion correction following Grimme */
1350 /* S.Grimme, J.Comput.Chem. 27, (2006), 1787-1799. */
1352 int
1353 DispersionCmd(ClientData clientData,
1354 Tcl_Interp *interp,
1355 int argc,
1356 const char *argv[])
1357 {
1359 /* Local variables */
1360 Frag coords;
1361 Matrix energy, gradient;
1362 ObjList lcoords, lene, lgrad;
1364 const char *s;
1365 int nat,ipar,iat,jat,iz,jz;
1367 /* string-buffers for use in Tcl_Eval , Tcl_SetVar
1368 (to avoid conflict with strbuff altered inside the called C-routines */
1369 char tclbuff[STRBUFFLEN];
1370 char tclbuff2[STRBUFFLEN];
1372 char buff2[2];
1373 #define MAXPAR 100
1374 double rnod[MAXPAR],c6[MAXPAR];
1375 double rsum,dist,c6ij,dexp;
1376 double edisp,fdamp,s6,grad,edamp,c_over_r6;
1377 struct atom_struct *a1, *a2;
1378 struct vector_struct v12, f12;
1379 int verbose; /* how much output to generate: 0 or 1 */
1380 int first; /* if 1 do not calculate, but test if parameters are set, write citation info */
1382 /* scaling factor for the funcational */
1383 GetVar(s,"s6");
1384 if(Tcl_GetDouble(interp,s,&s6)){
1385 printf("bad s6 in dispersion: %s\n",s);
1386 return TCL_ERROR;
1387 }
1389 /* verbosity */
1390 GetVar(s,"verbose");
1391 if(Tcl_GetInt(interp,s,&verbose)){
1392 printf("bad verbose in dispersion: %s\n",s);
1393 return TCL_ERROR;
1394 }
1396 GetVar(s,"first");
1397 if(Tcl_GetInt(interp,s,&first)){
1398 printf("bad first in dispersion: %s\n",s);
1399 return TCL_ERROR;
1400 }
1402 /* Access initial structure */
1403 GetVar(s,"coords");
1404 lcoords = get_objlist(s,"fragment",CHEMSH_OBJ_EXISTS,0);
1405 if(!lcoords){
1406 printf("bad coords tag in dispersion_correction: %s\n",s);
1407 return TCL_ERROR;
1408 }
1409 coords = (Frag) lcoords->object->data;
1410 nat = FRAG_get_number_of_atoms(coords);
1412 /* Gradient matrix to accumulate into (3n*1) */
1413 GetVar(s,"gradient");
1414 lgrad = get_objlist(s,"matrix",CHEMSH_OBJ_EXISTS,0);
1415 if(!lgrad){
1416 /* printf("bad matrix tag for gradient in dispersion_correction: %s \n",s);
1417 return TCL_ERROR; */
1418 gradient = NULL;
1419 } else {
1420 gradient = (Matrix) lgrad->object->data;
1421 }
1423 /* energy matrix (1*1) */
1424 GetVar(s,"energy");
1425 lene = get_objlist(s,"matrix",CHEMSH_OBJ_EXISTS,0);
1426 if(!lene){
1427 printf("bad matrix tag for energy in dispersion_correction: %s\n", s);
1428 return TCL_ERROR;
1429 }
1430 energy = (Matrix) lene->object->data;
1432 /* set parameters according to Table 1 in
1433 S.Grimme, J.Comput.Chem. 27, (2006), 1787-1799.
1435 rnod is R0, the vdW radius
1436 c6 is the C6 coefficient per atom
1437 rsum is the summ of the vdW radii of the two atoms
1438 dist is the distance between the atoms
1439 c6ij is the geometric mean of the c6 parameters of the two atoms
1440 dexp is the exponential factor d in Eq. 12 in Grimme's paper
1441 */
1443 /* set default values for atom types that are not parametrised */
1444 for(ipar=1;ipar<MAXPAR;ipar++){
1445 rnod[ipar]=1.e+99;
1446 c6[ipar]=0.;
1447 }
1449 dexp=20.;
1451 if(first){
1452 printf("\nAdding a dispersion correction to the energy. Please cite:\n");
1453 printf("S. Grimme, J. Comput. Chem. 27, (2006), 1787-1799\n\n");
1454 printf("Functional-dependent prefactor (s6): %f\n",s6);
1455 }
1457 rnod[1]= 1.001 ;
1458 rnod[2]= 1.012 ;
1459 rnod[3]= 0.825 ;
1460 rnod[4]= 1.408 ;
1461 rnod[5]= 1.485 ;
1462 rnod[6]= 1.452 ;
1463 rnod[7]= 1.397 ;
1464 rnod[8]= 1.342 ;
1465 rnod[9]= 1.287 ;
1466 rnod[10]= 1.243 ;
1467 rnod[11]= 1.144 ;
1468 rnod[12]= 1.364 ;
1469 rnod[13]= 1.639 ;
1470 rnod[14]= 1.716 ;
1471 rnod[15]= 1.705 ;
1472 rnod[16]= 1.683 ;
1473 rnod[17]= 1.639 ;
1474 rnod[18]= 1.595 ;
1475 rnod[19]= 1.485 ;
1476 rnod[20]= 1.474 ;
1477 rnod[21]= 1.562 ;
1478 rnod[22]= 1.562 ;
1479 rnod[23]= 1.562 ;
1480 rnod[24]= 1.562 ;
1481 rnod[25]= 1.562 ;
1482 rnod[26]= 1.562 ;
1483 rnod[27]= 1.562 ;
1484 rnod[28]= 1.562 ;
1485 rnod[29]= 1.562 ;
1486 rnod[30]= 1.562 ;
1487 rnod[31]= 1.650 ;
1488 rnod[32]= 1.727 ;
1489 rnod[33]= 1.760 ;
1490 rnod[34]= 1.771 ;
1491 rnod[35]= 1.749 ;
1492 rnod[36]= 1.727 ;
1493 rnod[37]= 1.628 ;
1494 rnod[38]= 1.606 ;
1495 rnod[39]= 1.639 ;
1496 rnod[40]= 1.639 ;
1497 rnod[41]= 1.639 ;
1498 rnod[42]= 1.639 ;
1499 rnod[43]= 1.639 ;
1500 rnod[44]= 1.639 ;
1501 rnod[45]= 1.639 ;
1502 rnod[46]= 1.639 ;
1503 rnod[47]= 1.639 ;
1504 rnod[48]= 1.639 ;
1505 rnod[49]= 1.672 ;
1506 rnod[50]= 1.804 ;
1507 rnod[51]= 1.881 ;
1508 rnod[52]= 1.892 ;
1509 rnod[53]= 1.892 ;
1510 rnod[54]= 1.881 ;
1512 c6[1]= 0.14 ;
1513 c6[2]= 0.08 ;
1514 c6[3]= 1.61 ;
1515 c6[4]= 1.61 ;
1516 c6[5]= 3.13 ;
1517 c6[6]= 1.75 ;
1518 c6[7]= 1.23 ;
1519 c6[8]= 0.70 ;
1520 c6[9]= 0.75 ;
1521 c6[10]= 0.63 ;
1522 c6[11]= 5.71 ;
1523 c6[12]= 5.71 ;
1524 c6[13]= 10.79 ;
1525 c6[14]= 9.23 ;
1526 c6[15]= 7.84 ;
1527 c6[16]= 5.57 ;
1528 c6[17]= 5.07 ;
1529 c6[18]= 4.61 ;
1530 c6[19]= 10.80 ;
1531 c6[20]= 10.80 ;
1532 c6[21]= 10.80 ;
1533 c6[22]= 10.80 ;
1534 c6[23]= 10.80 ;
1535 c6[24]= 10.80 ;
1536 c6[25]= 10.80 ;
1537 c6[26]= 10.80 ;
1538 c6[27]= 10.80 ;
1539 c6[28]= 10.80 ;
1540 c6[29]= 10.80 ;
1541 c6[30]= 10.80 ;
1542 c6[31]= 16.99 ;
1543 c6[32]= 17.10 ;
1544 c6[33]= 16.37 ;
1545 c6[34]= 12.64 ;
1546 c6[35]= 12.47 ;
1547 c6[36]= 12.01 ;
1548 c6[37]= 24.67 ;
1549 c6[38]= 24.67 ;
1550 c6[39]= 24.67 ;
1551 c6[40]= 24.67 ;
1552 c6[41]= 24.67 ;
1553 c6[42]= 24.67 ;
1554 c6[43]= 24.67 ;
1555 c6[44]= 24.67 ;
1556 c6[45]= 24.67 ;
1557 c6[46]= 24.67 ;
1558 c6[47]= 24.67 ;
1559 c6[48]= 24.67 ;
1560 c6[49]= 37.32 ;
1561 c6[50]= 38.71 ;
1562 c6[51]= 38.44 ;
1563 c6[52]= 31.74 ;
1564 c6[53]= 31.50 ;
1565 c6[54]= 29.99 ;
1567 edisp=0.;
1569 for(iat=0;iat<nat;iat++){
1570 a1 = &coords->atoms[iat];
1571 iz = a1->znum;
1572 if(iz==0){
1573 printf("Error: could not determine Z of atom %d\n", iat+1);
1574 return TCL_ERROR;
1575 }
1576 if(iz>100){
1577 printf("Error: dispersion correction for atom %d (Element %d) not possible!\n",
1578 iat+1,iz);
1579 return TCL_ERROR;
1580 }
1581 if(iz>54 && first){
1582 printf("Warning: dispersion correction for atom %d (Element %d) not parametrised!\n",
1583 iat+1,iz);
1584 }
1585 for(jat=iat+1;jat<nat;jat++){
1587 a2 = &coords->atoms[jat];
1588 jz= a2->znum;
1590 if(jz==0){
1591 printf("Error: j could not determine Z of atom %d\n", jat+1);
1592 return TCL_ERROR;
1593 }
1594 if(jz>100){
1595 printf("Error: dispersion correction for atom %d (Element %d) not parametrised!\n",
1596 jat+1,jz);
1597 return TCL_ERROR;
1598 }
1599 if(jz>54 && first){
1600 printf("Warning: dispersion correction for atom %d (Element %d) not implemented!\n",
1601 jat+1,jz);
1602 }
1604 /* Implementation of Equations 11, 12, and 13 of the Grimme paper */
1607 rsum = (rnod[iz] + rnod[jz] ) / 0.52917706 ; /* convert ang to au */
1608 c6ij=sqrt( c6[iz] * c6[jz] ) * 17.345304 ; /* convert J*nm^6/mol to au */
1610 a1 = &coords->atoms[iat];
1611 a2 = &coords->atoms[jat];
1612 dist = vdist(a1->pos,a2->pos);
1614 v12 = vnorm(vdiff(a2->pos,a1->pos));
1616 edamp = exp( -dexp * (dist/rsum - 1.) ) ;
1618 c_over_r6 = c6ij / ( dist *dist *dist *dist *dist *dist);
1620 fdamp=1./( 1. + edamp );
1622 edisp += -s6 * c_over_r6 * fdamp;
1624 if(gradient) {
1625 if(gradient->d && !first){
1626 grad = s6 * ( -6. * c_over_r6 * fdamp / dist +
1627 c_over_r6 * dexp * edamp / rsum /
1628 ( 1. + edamp ) / ( 1. + edamp ) );
1629 f12 = vscal(v12,grad);
1630 gradient->d[iat*3 + 0] += f12.x[0];
1631 gradient->d[iat*3 + 1] += f12.x[1];
1632 gradient->d[iat*3 + 2] += f12.x[2];
1633 gradient->d[jat*3 + 0] -= f12.x[0];
1634 gradient->d[jat*3 + 1] -= f12.x[1];
1635 gradient->d[jat*3 + 2] -= f12.x[2];
1636 }
1637 }
1639 }
1640 }
1643 if( ! first ) {
1645 energy->d[0] += edisp ;
1647 if(verbose >= 1) printf("Dispersion correction included in the energy: %10.5f Hartree\n",edisp);
1649 }
1651 rel_objlist(lcoords);
1652 rel_objlist(lene);
1653 if(gradient) rel_objlist(lgrad);
1656 return TCL_OK;
1658 }
