/****************************************************************************** COPYRIGHT - Copyright and any other intellectual property rights in this software belong to CLRC. the recipient may not copy the software except for backup purposes. Nor may the recipient re-distribute the software to any third party without the prior written consent of CLRC. WARRANTIES - This software is supplied on an "AS IS" basis and no representations are made, nor warranties are given as to its accuracy, usefulness or performance. Nor do CLRC give any warranties, express or implied, as to its merchantability or fitness for a particular purpose. LIABILITY - CLRC will not be liable for an loss of any kind arising from the supply and use of this software. **************************************************************************** $Author: tbenighaus $ $Date: 2010-09-23 09:35:53 +0200 (Thu, 23. Sep 2010) $ $Revision: 3689 $ $Source$ $State$ **************************************************************************** C code for the restraints function Providing a flexible way to augment any energy/gradient evaluator with some additional terms, typically for umbrella sampling and related methods */ #include #include #include #include #include #include #include "tcl.h" #include "objects.h" #include "chemsh.h" #include "dutil.h" #include "dblock.h" #include "dmatrix.h" #include "dfragment.h" #include "restraint.h" #if defined (_WIN32) || defined (CRAYXX) #define CAPITALSYMBOLS #endif #ifdef CAPITALSYMBOLS #define torsionval_ TORSIONVAL #define torsionderivative_ TORSIONDERIVATIVE #endif #ifdef NOC_ #define torsionval_ torsionval #define torsionderivative_ torsionderivative #endif extern Tcl_Interp *chemsh_interp; static int debug = 0; /* Declarations */ int rst_newconfig(Tcl_Interp *interp,struct rst_struct *rst, int argc,const char **argv); int calc_rst_ff(Ff ff, Frag f, Matrix e, Matrix g); int rst_set_config(Tcl_Interp *interp, struct rst_struct *rst); Tcl_CmdProc Chemsh_RestraintObjCmd, RestraintCmd, DispersionCmd; void torsionval_(double torpos[4][3], double *phi); void torsionderivative_(double torpos[4][3], double *phi, double torder[4][3]); void declare_restraint(Tcl_Interp *interp){ Tcl_CreateCommand(interp, "restraint_c" ,RestraintCmd, (ClientData) NULL,(Tcl_CmdDeleteProc *) NULL); Tcl_CreateCommand(interp, "dispersion_c" ,DispersionCmd, (ClientData) NULL,(Tcl_CmdDeleteProc *) NULL); } int RestraintCmd(ClientData clientData, Tcl_Interp *interp, int argc, const char *argv[]) { int Chemsh_RestraintObjCmd(); /* Local variables */ struct rst_struct *rst; const char *s, **work, **work2; int itmp, i, j; double atmp; const char *cmd; int nfld; Ff ff; int ff_counter; /* string-buffers for use in Tcl_Eval , Tcl_SetVar (to avoid conflict with strbuff altered inside the called C-routines */ char tclbuff[STRBUFFLEN]; char tclbuff2[STRBUFFLEN]; char buff2[2]; cmd = argv[0]; rst = (struct rst_struct *)ckalloc(sizeof(struct rst_struct)); rst->current_ff = NULL; /* store object name */ strcpym(&rst->name, argv[1]); /* create the Tcl interface to the command */ if(debug)printf("create object %s\n",argv[1]); Tcl_CreateCommand(interp, rst->name, Chemsh_RestraintObjCmd, (ClientData) rst, (Tcl_CmdDeleteProc *) NULL); /* Access initial structure */ GetArg(s,"coords","restraints"); rst->lc = get_objlist(s,"fragment",CHEMSH_OBJ_EXISTS,0); if(!rst->lc){ printf("bad tag coords=%s\n",s); return TCL_ERROR; } rst->c = (Frag) rst->lc->object->data; rst->natms = FRAG_get_number_of_atoms(rst->c); /* Gradient matrix to accumulate into (3n*1) */ GetArg(s,"gradient","restraints"); rst->lg = get_objlist(s,"matrix",CHEMSH_OBJ_EXISTS,0); if(!rst->lg){ printf("%s: bad matrix tag %s\n", cmd,s); return TCL_ERROR; } rst->g = (Matrix) rst->lg->object->data; /* energy matrix (1*1) */ GetArg(s,"energy","restraints"); rst->le = get_objlist(s,"matrix",CHEMSH_OBJ_EXISTS,0); if(!rst->le){ printf("%s: bad matrix tag %s\n", cmd,s); return TCL_ERROR; } rst->e = (Matrix) rst->le->object->data; /* Now parse requests for terms Atoms are identified by integer labels (counting from 1) */ sprintf(tclbuff,"rst_parm(%s,%s)",rst->name,"restraints"); if((s = Tcl_GetVar(interp,tclbuff,0)) == NULL){ printf("no restraints setting in restraints module\n"); return TCL_ERROR; } if(debug)printf("restraints: %s\n",s); Tcl_SplitList(interp,s,&rst->nrest,&work); if(debug)printf("parsing %d restraints\n",rst->nrest); for(j=0;jnrest;j++){ Tcl_SplitList(interp,work[j],&nfld,&work2); ff=new_rst_ff(rst); ff->index = j+1; if(!strcmp(work2[0],"bond")){ /* harmonic distance restraint - expected format is bond i j r0 K */ if ( nfld != 5 ) { printf("restraint: invalid format for restraint bond\n"); printf("number of field exprected: 5 (bond atom1 atom2 r0 K) \n"); printf("number of fields found: %d\n",nfld); return TCL_ERROR; } if(Tcl_GetInt(interp,work2[1],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]); return TCL_ERROR; } ff->i1 = itmp-1; if(Tcl_GetInt(interp,work2[2],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]); return TCL_ERROR; } ff->i2 = itmp-1; ff->type = BOND; if(Tcl_GetDouble(interp,work2[3],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]); return TCL_ERROR; } ff->r0 = atmp; if(Tcl_GetDouble(interp,work2[4],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]); return TCL_ERROR; } ff->fc = atmp; } else if(!strcmp(work2[0],"elec")){ /* electrostatic term - expected format is elec i j chargeprod not a restraint, but added in as it may be useful for other purposes..... */ if ( nfld != 4 ) { printf("restraint: invalid format for restraint elec\n"); printf("number of field exprected: 4 (elec i j chargeprod) \n"); printf("number of fields found: %d\n",nfld); return TCL_ERROR; } if(Tcl_GetInt(interp,work2[1],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]); return TCL_ERROR; } ff->i1 = itmp-1; if(Tcl_GetInt(interp,work2[2],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]); return TCL_ERROR; } ff->i2 = itmp-1; ff->type = ELECTROSTATIC; if(Tcl_GetDouble(interp,work2[3],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]); return TCL_ERROR; } ff->chargeprod = atmp; }else if(!strcmp(work2[0],"angle")){ /* harmonic angle restraint - expected format is angle i j k r0 K */ if ( nfld != 6 ) { printf("restraint: invalid format for restraint angle\n"); printf("number of field exprected: 6 (angle i j k r0 K) \n"); printf("number of fields found: %d\n",nfld); return TCL_ERROR; } if(Tcl_GetInt(interp,work2[1],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]); return TCL_ERROR; } ff->i1 = itmp-1; if(Tcl_GetInt(interp,work2[2],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]); return TCL_ERROR; } ff->i2 = itmp-1; if(Tcl_GetInt(interp,work2[3],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]); return TCL_ERROR; } ff->i3 = itmp-1; ff->type = ANG; if(Tcl_GetDouble(interp,work2[4],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]); return TCL_ERROR; } /* Store internally r0 in radians */ #define cnv 57.29577951 ff->r0 = atmp / cnv; if(Tcl_GetDouble(interp,work2[5],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]); return TCL_ERROR; } ff->fc = atmp; }else if(!strcmp(work2[0],"torsion")){ /* harmonic torsion restraint - expected format is torsion i j k l r0 K */ if ( nfld != 7 ) { printf("restraint: invalid format for restraint torsion\n"); printf("number of field exprected: 7 (torsion i j k l r0 K) \n"); printf("number of fields found: %d\n",nfld); return TCL_ERROR; } if(Tcl_GetInt(interp,work2[1],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]); return TCL_ERROR; } ff->i1 = itmp-1; if(Tcl_GetInt(interp,work2[2],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]); return TCL_ERROR; } ff->i2 = itmp-1; if(Tcl_GetInt(interp,work2[3],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]); return TCL_ERROR; } ff->i3 = itmp-1; if(Tcl_GetInt(interp,work2[4],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]); return TCL_ERROR; } ff->i4 = itmp-1; ff->type = TOR; if(Tcl_GetDouble(interp,work2[5],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]); return TCL_ERROR; } /* Store internally r0 in radians */ ff->r0 = atmp / cnv; if(Tcl_GetDouble(interp,work2[6],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",6,work2[6]); return TCL_ERROR; } ff->fc = atmp; }else if(!strcmp(work2[0],"bonddiff2")){ /* difference between two bond lengths restraint - expected format is bonddiff2 i j k l r0 K */ if ( nfld != 7 ) { printf("restraint: invalid format for restraint bonddiff2\n"); printf("number of field exprected: 7 (bonddiff2 i j k l r0 K) \n"); printf("number of fields found: %d\n",nfld); return TCL_ERROR; } if(Tcl_GetInt(interp,work2[1],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]); return TCL_ERROR; } ff->i1 = itmp-1; if(Tcl_GetInt(interp,work2[2],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]); return TCL_ERROR; } ff->i2 = itmp-1; if(Tcl_GetInt(interp,work2[3],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]); return TCL_ERROR; } ff->i3 = itmp-1; if(Tcl_GetInt(interp,work2[4],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]); return TCL_ERROR; } ff->i4 = itmp-1; ff->type = BONDDIFF2; if(Tcl_GetDouble(interp,work2[5],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]); return TCL_ERROR; } ff->r0 = atmp; if(Tcl_GetDouble(interp,work2[6],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",6,work2[6]); return TCL_ERROR; } ff->fc = atmp; }else if(!strcmp(work2[0],"bonddiff3")){ /* difference between three bond lengths restraint - expected format is bonddiff3 i j k l m n r0 K */ if ( nfld != 9 ) { printf("restraint: invalid format for restraint bonddiff3\n"); printf("number of field exprected: 9 (bonddiff3 i j k l m n r0 K) \n"); printf("number of fields found: %d\n",nfld); return TCL_ERROR; } if(Tcl_GetInt(interp,work2[1],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]); return TCL_ERROR; } ff->i1 = itmp-1; if(Tcl_GetInt(interp,work2[2],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]); return TCL_ERROR; } ff->i2 = itmp-1; if(Tcl_GetInt(interp,work2[3],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]); return TCL_ERROR; } ff->i3 = itmp-1; if(Tcl_GetInt(interp,work2[4],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]); return TCL_ERROR; } ff->i4 = itmp-1; if(Tcl_GetInt(interp,work2[5],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]); return TCL_ERROR; } ff->i5 = itmp-1; if(Tcl_GetInt(interp,work2[6],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",6,work2[6]); return TCL_ERROR; } ff->i6 = itmp-1; ff->type = BONDDIFF3; if(Tcl_GetDouble(interp,work2[7],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",7,work2[7]); return TCL_ERROR; } ff->r0 = atmp; if(Tcl_GetDouble(interp,work2[8],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",8,work2[8]); return TCL_ERROR; } ff->fc = atmp; }else if(!strcmp(work2[0],"bonddiff4")){ /* difference between four bond lengths restraint - expected format is bonddiff3 i j k l m n o p r0 K E=K/2 * ( d(ij)+d(kl)-d(mn)-d(op) - r0 )**2 */ if ( nfld != 11 ) { printf("restraint: invalid format for restraint bonddiff4\n"); printf("number of field exprected: 11 (bonddiff3 i j k l m n o p r0 K) \n"); printf("number of fields found: %d\n",nfld); return TCL_ERROR; } if(Tcl_GetInt(interp,work2[1],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]); return TCL_ERROR; } ff->i1 = itmp-1; if(Tcl_GetInt(interp,work2[2],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]); return TCL_ERROR; } ff->i2 = itmp-1; if(Tcl_GetInt(interp,work2[3],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]); return TCL_ERROR; } ff->i3 = itmp-1; if(Tcl_GetInt(interp,work2[4],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]); return TCL_ERROR; } ff->i4 = itmp-1; if(Tcl_GetInt(interp,work2[5],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]); return TCL_ERROR; } ff->i5 = itmp-1; if(Tcl_GetInt(interp,work2[6],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",6,work2[6]); return TCL_ERROR; } ff->i6 = itmp-1; if(Tcl_GetInt(interp,work2[7],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",7,work2[7]); return TCL_ERROR; } ff->i7 = itmp-1; if(Tcl_GetInt(interp,work2[8],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",8,work2[8]); return TCL_ERROR; } ff->i8 = itmp-1; ff->type = BONDDIFF4; if(Tcl_GetDouble(interp,work2[9],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",7,work2[9]); return TCL_ERROR; } ff->r0 = atmp; if(Tcl_GetDouble(interp,work2[10],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",8,work2[10]); return TCL_ERROR; } ff->fc = atmp; }else if(!strcmp(work2[0],"sbc_atom")){ /* Spherical boundary conditions around an atom - expected format is sbc_atom i r0 K E= Sum ( 0.5 * K * (| r_a - r_i | - r0 )^2 ) */ if ( nfld != 4 ) { printf("restraint: invalid format for restraint sbc_atom\n"); printf("number of field exprected: 4 (sbc_atom i r0 K) \n"); printf("number of fields found: %d\n",nfld); return TCL_ERROR; } if(Tcl_GetInt(interp,work2[1],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]); return TCL_ERROR; } ff->i1 = itmp-1; if(Tcl_GetDouble(interp,work2[2],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]); return TCL_ERROR; } ff->r0 = atmp; if(Tcl_GetDouble(interp,work2[3],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]); return TCL_ERROR; } ff->fc = atmp; ff->type = SBC_ATOM; }else if(!strcmp(work2[0],"sbc")){ /* Spherical boundary conditions around an geometric center c, specified by x,y,z - expected format is sbc x y z r0 K E= Sum ( 0.5 * K * (| r_a - c | - r0 )^2 ) */ if ( nfld != 6 ) { printf("restraint: invalid format for restraint sbc_atom\n"); printf("number of field exprected: 6 (sbc x y z r0 K) \n"); printf("number of fields found: %d\n",nfld); return TCL_ERROR; } if(Tcl_GetDouble(interp,work2[1],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]); return TCL_ERROR; } ff->x = atmp; if(Tcl_GetDouble(interp,work2[2],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]); return TCL_ERROR; } ff->y = atmp; if(Tcl_GetDouble(interp,work2[3],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]); return TCL_ERROR; } ff->z = atmp; if(Tcl_GetDouble(interp,work2[4],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]); return TCL_ERROR; } ff->r0 = atmp; if(Tcl_GetDouble(interp,work2[5],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]); return TCL_ERROR; } ff->fc = atmp; ff->type = SBC; }else if(!strcmp(work2[0],"sbc_individual")){ /* Spherical boundary conditions around an geometric center c for an explicit atom, specified by x,y,z - expected format is sbc_individual atom_number x y z r0 K E= Sum ( 0.5 * K * (| r_a - c | - r0 )^2 ) */ if ( nfld != 7 ) { printf("restraint: invalid format for restraint sbc_individual\n"); printf("number of field exprected: 7 (sbc_individual atom_number x y z r0 K) \n"); printf("number of fields found: %d\n",nfld); return TCL_ERROR; } if(Tcl_GetInt(interp,work2[1],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]); return TCL_ERROR; } ff->i1 = itmp-1; if(Tcl_GetDouble(interp,work2[2],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]); return TCL_ERROR; } ff->x = atmp; if(Tcl_GetDouble(interp,work2[3],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]); return TCL_ERROR; } ff->y = atmp; if(Tcl_GetDouble(interp,work2[4],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]); return TCL_ERROR; } ff->z = atmp; if(Tcl_GetDouble(interp,work2[5],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]); return TCL_ERROR; } ff->r0 = atmp; if(Tcl_GetDouble(interp,work2[6],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]); return TCL_ERROR; } ff->fc = atmp; ff->type = SBC_IND; }else if(!strcmp(work2[0],"cart_atom")){ /* Restrains an atom around its original position cart_atom iat K E= 0.5 * K * (| r_a - c | )^2 */ if ( nfld != 3 ) { printf("restraint: invalid format for restraint cart_atom\n"); printf("number of field exprected: 3 (cart_atom atom_number K) \n"); printf("number of fields found: %d\n",nfld); return TCL_ERROR; } if(Tcl_GetInt(interp,work2[1],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]); return TCL_ERROR; } ff->i1 = itmp-1; ff->x = rst->c->atoms[itmp-1].pos.x[0]; ff->y = rst->c->atoms[itmp-1].pos.x[1]; ff->z = rst->c->atoms[itmp-1].pos.x[2]; ff->r0 = 0.; if(Tcl_GetDouble(interp,work2[2],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]); return TCL_ERROR; } ff->fc = atmp; ff->type = CART_ATOM; }else if(!strcmp(work2[0],"cart")){ /* Restrains an atom around the position x y z cart iat x y z K E= 0.5 * K * (| r_a - c | )^2 */ if ( nfld != 6 ) { printf("restraint: invalid format for restraint cart\n"); printf("number of field exprected: 6 (cart atom_number x y z K) \n"); printf("number of fields found: %d\n",nfld); return TCL_ERROR; } if(Tcl_GetInt(interp,work2[1],&itmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",1,work2[1]); return TCL_ERROR; } ff->i1 = itmp-1; if(Tcl_GetDouble(interp,work2[2],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",2,work2[2]); return TCL_ERROR; } ff->x = atmp; if(Tcl_GetDouble(interp,work2[3],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",3,work2[3]); return TCL_ERROR; } ff->y = atmp; if(Tcl_GetDouble(interp,work2[4],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",4,work2[4]); return TCL_ERROR; } ff->z = atmp; ff->r0 = 0.; if(Tcl_GetDouble(interp,work2[5],&atmp)){ printf("restraint: invalid value for restraint parameter field %d (%s)\n",5,work2[5]); return TCL_ERROR; } ff->fc = atmp; ff->type = CART; } else { printf("restraint: invalid value for restraint type (%s)\n",work2[0]); return TCL_ERROR; } ckfree((char *) work2); } ckfree((char *) work); printf("Restraints Atoms (object name %s) :\n",rst->name); for(ff=rst->last_ff;ff;ff=ff->prev){ switch(ff->type){ case BOND: printf("Harmonic Distance %5d %5d r0 = %f fc = %f \n", ff->i1+1, ff->i2+1, ff->r0, ff->fc ); break; case ANG: printf("Harmonic Angle %5d %5d %5d r0 = %f fc = %f \n", ff->i1+1, ff->i2+1, ff->i3+1, ff->r0, ff->fc ); break; case ELECTROSTATIC: printf("Electrostatic %5d %5d qiqj = %f \n", ff->i1+1, ff->i2+1, ff->chargeprod); break; case TOR: printf("Torsion %5d %5d %5d %5d r0 = %f fc = %f \n", ff->i1+1, ff->i2+1, ff->i3+1, ff->i4+1, ff->r0, ff->fc ); break; case BONDDIFF2: printf("Bonddiff-2 %5d %5d %5d %5d r0 = %f fc = %f \n", ff->i1+1, ff->i2+1, ff->i3+1, ff->i4+1, ff->r0, ff->fc ); break; case BONDDIFF3: printf("Bonddiff-3 %5d %5d %5d %5d %5d %5d r0 = %f fc = %f \n", ff->i1+1, ff->i2+1, ff->i3+1, ff->i4+1, ff->i5+1, ff->i6+1, ff->r0, ff->fc ); break; case BONDDIFF4: printf("Bonddiff-4 %5d %5d %5d %5d %5d %5d %5d %5d r0 = %f fc = %f \n", ff->i1+1, ff->i2+1, ff->i3+1, ff->i4+1, ff->i5+1, ff->i6+1, ff->i7+1, ff->i8+1, ff->r0, ff->fc ); break; case SBC_ATOM: printf("SBC Atom %5d r0 = %f fc = %f \n", ff->i1+1, ff->r0, ff->fc ); break; case SBC: printf("SBC %f %f %f r0 = %f fc = %f \n", ff->x, ff->y, ff->z, ff->r0, ff->fc ); break; case SBC_IND: printf("SBC_individual %5d %f %f %f r0 = %f fc = %f \n", ff->i1+1, ff->x, ff->y, ff->z, ff->r0, ff->fc ); break; case CART_ATOM: printf("Cart Atom %5d %10.5f %10.5f %10.5f fc = %f \n", ff->i1+1, ff->x, ff->y, ff->z, ff->fc ); break; case CART: printf("Cart %5d %10.5f %10.5f %10.5f fc = %f \n", ff->i1+1, ff->x, ff->y, ff->z, ff->fc ); break; } } return TCL_OK; } /* Implements the commands associated with an restraint object The first argument is the required method */ int Chemsh_RestraintObjCmd(ClientData clientData, Tcl_Interp *interp, int argc, const char *argv[]) { /* local variables */ char tclbuff[STRBUFFLEN]; char strerr[5]; char c; int i, j, length; struct rst_struct *rst; Ff ff; int itmp; rst = (struct rst_struct *) clientData; /* default Tcl result */ strcpy(strerr,"0"); if (argc == 1){ printf("optimiser object %s invoked without a command\n",rst->name); return TCL_ERROR; } c = argv[1][0]; /* First letter of command */ length = strlen(argv[1]); sprintf(tclbuff,"global rst_parm"); if (Tcl_Eval(interp,tclbuff) != TCL_OK) return TCL_ERROR; if(debug) printf("%s Invoking method: %s\n",rst->name, argv[1]); if ((c == 'd') && ( strncmp(argv[1],"destroy", length) == 0)) { /* nothing to do really... maybe release rst structure ??? */ rel_objlist(rst->lc); rel_objlist(rst->lg); rel_objlist(rst->le); /* delete_restraint_memory(rst); */ /*pop_module_name(interp); */ } else if ((c == 'c') && ( strncmp(argv[1],"configure", length) == 0)) { /* Configure */ if(rst_newconfig(interp,rst,argc-2,argv+2)) { return TCL_ERROR; } /*rst_print_parms(rst);*/ } else if ((c == 'o') && ( strncmp(argv[1],"output", length) == 0)) { /* Print current values of restraint values and energies and force terms */ } else if ((c == 'a') && ( strncmp(argv[1],"accumulate", length) == 0)) { /* Accumulate and print results */ for(ff=rst->last_ff;ff;ff=ff->prev){ calc_rst_ff(ff, rst->c, rst->e, rst->g); } } else { printf("unrecognised method for restraint object %s : %s\n",rst->name, argv[1]); sprintf(tclbuff,"unrecognised method for restraint object %s : %s\n",rst->name, argv[1]); Tcl_SetResult(interp, tclbuff, TCL_VOLATILE); return TCL_ERROR; } /* clean exits */ Tcl_SetResult(interp, strerr, TCL_VOLATILE); return TCL_OK; } /* Process configuration options */ int rst_newconfig(Tcl_Interp *interp, struct rst_struct *rst, int argc, const char **argv) { char tclbuff[STRBUFFLEN]; int i; for(i=0;iname); if(Tcl_Eval(interp,tclbuff)){ printf("%s: configure failed - %s\n",rst->name,interp->result); return TCL_ERROR; } rst_set_config(interp,rst); printf("%s: configure completed: ",rst->name); for(i=0;itype){ case BOND: a1 = &f->atoms[ff->i1]; a2 = &f->atoms[ff->i2]; r12 = vdist(a1->pos,a2->pos); v12 = vnorm(vdiff(a2->pos,a1->pos)); d = r12 - ff->r0; f12 = vscal(v12,-d*ff->fc); if(g->d){ g->d[ff->i1*3 + 0] += f12.x[0]; g->d[ff->i1*3 + 1] += f12.x[1]; g->d[ff->i1*3 + 2] += f12.x[2]; g->d[ff->i2*3 + 0] -= f12.x[0]; g->d[ff->i2*3 + 1] -= f12.x[1]; g->d[ff->i2*3 + 2] -= f12.x[2]; } ee = 0.5*d*d*ff->fc; e->d[0] += ee; if(debug)printf("r,energy in calc_ff %lf %10.4e\n",r12,ee); printf("Restraint %2d dist : value %11.4e energy %11.4e force %11.4e\n",ff->index,r12,ee,-d*ff->fc); break; case ELECTROSTATIC: a1 = &f->atoms[ff->i1]; a2 = &f->atoms[ff->i2]; r12 = vdist(a1->pos,a2->pos); v12 = vnorm(vdiff(a2->pos,a1->pos)); f12 = vscal(v12,ff->chargeprod*(1.0/(r12*r12))); if(g->d){ g->d[ff->i1*3 + 0] += f12.x[0]; g->d[ff->i1*3 + 1] += f12.x[1]; g->d[ff->i1*3 + 2] += f12.x[2]; g->d[ff->i2*3 + 0] -= f12.x[0]; g->d[ff->i2*3 + 1] -= f12.x[1]; g->d[ff->i2*3 + 2] -= f12.x[2]; } ee = ff->chargeprod*(1.0/r12); e->d[0] += ee; if(debug)printf("chargeprod, r,energy in calc_ff %f %lf %10.4e\n",ff->chargeprod,r12,ee); printf("Restraint %2d elec : value %11.4e energy %11.4e\n",ff->index,r12,ee); break; case ANG: a1 = &f->atoms[ff->i1]; a2 = &f->atoms[ff->i2]; a3 = &f->atoms[ff->i3]; r12 = vdist(a1->pos,a2->pos); r32 = vdist(a3->pos,a2->pos); v12 = vnorm(vdiff(a2->pos,a1->pos)); v32 = vnorm(vdiff(a2->pos,a3->pos)); cosphi = vdot(v12,v32); phi = acos(cosphi); sinphi = sin(phi); d = phi - ff->r0; f1 = vscal(vdiff(vscal(v12,cosphi),v32), d*ff->fc* 1.0/(r12*sinphi)); f3 = vscal(vdiff(vscal(v32,cosphi),v12), d*ff->fc* 1.0/(r32*sinphi)); if(debug){ vprint("v12",v12); vprint("v32",v32); vprint("f1",f1); vprint("f3",f3); } if(g->d){ g->d[ff->i1*3 + 0] -= f1.x[0]; g->d[ff->i1*3 + 1] -= f1.x[1]; g->d[ff->i1*3 + 2] -= f1.x[2]; g->d[ff->i3*3 + 0] -= f3.x[0]; g->d[ff->i3*3 + 1] -= f3.x[1]; g->d[ff->i3*3 + 2] -= f3.x[2]; g->d[ff->i2*3 + 0] += (f1.x[0] + f3.x[0]); g->d[ff->i2*3 + 1] += (f1.x[1] + f3.x[1]); g->d[ff->i2*3 + 2] += (f1.x[2] + f3.x[2]); } ee = 0.5 * d * d * ff->fc; if(debug)printf("phi,energy in calc_ff %lf %10.4e\n",phi,ee); printf("Restraint %2d angl : value %11.4e energy %11.4e force %11.4e\n",ff->index,phi,ee,-d*ff->fc); e->d[0] += ee; break; case TOR: /* map to fortran routines in the dynamic-module */ for(icor=1;icor<=3;icor++) { a1 = &f->atoms[ff->i1]; torpos[0][icor-1]=a1->pos.x[icor-1]; a1 = &f->atoms[ff->i2]; torpos[1][icor-1]=a1->pos.x[icor-1]; a1 = &f->atoms[ff->i3]; torpos[2][icor-1]=a1->pos.x[icor-1]; a1 = &f->atoms[ff->i4]; torpos[3][icor-1]=a1->pos.x[icor-1]; } if(g->d){ torsionderivative_(torpos,&phi,torder); } else { torsionval_(torpos,&phi); } d = phi - ff->r0; /* check for a 360 deg rotation and avoid it */ pi = 4. * atan(1.); if( d > 7. || d < -7.) { printf("Error: torsion values differ too much!\n"); printf("Current value %f, centre of the restraint %f\n",phi,ff->r0); exit(1); } if(d > 3.1415 ) d=d-2.*pi; if(d < -3.1415 ) d=d+2.*pi; if(g->d){ fact= ff->fc * d; g->d[ff->i1*3 + 0] += fact*torder[0][0]; g->d[ff->i1*3 + 1] += fact*torder[0][1]; g->d[ff->i1*3 + 2] += fact*torder[0][2]; g->d[ff->i2*3 + 0] += fact*torder[1][0]; g->d[ff->i2*3 + 1] += fact*torder[1][1]; g->d[ff->i2*3 + 2] += fact*torder[1][2]; g->d[ff->i3*3 + 0] += fact*torder[2][0]; g->d[ff->i3*3 + 1] += fact*torder[2][1]; g->d[ff->i3*3 + 2] += fact*torder[2][2]; g->d[ff->i4*3 + 0] += fact*torder[3][0]; g->d[ff->i4*3 + 1] += fact*torder[3][1]; g->d[ff->i4*3 + 2] += fact*torder[3][2]; } ee = 0.5 * d * d * ff->fc; if(debug)printf("Torsion: phi, energy in calc_ff %lf %10.4e\n",phi,ee); printf("Restraint %2d tors : value %11.4e energy %11.4e force %11.4e\n",ff->index,phi,ee,-d*ff->fc); e->d[0] += ee; break; case BONDDIFF2: a1 = &f->atoms[ff->i1]; a2 = &f->atoms[ff->i2]; a3 = &f->atoms[ff->i3]; a4 = &f->atoms[ff->i4]; r12 = vdist(a1->pos,a2->pos); r34 = vdist(a3->pos,a4->pos); v12 = vnorm(vdiff(a2->pos,a1->pos)); v34 = vnorm(vdiff(a4->pos,a3->pos)); d = r12 - r34 - ff->r0; f12 = vscal(v12,-d*ff->fc); f34 = vscal(v34,-d*ff->fc); if(g->d){ g->d[ff->i1*3 + 0] += f12.x[0]; g->d[ff->i1*3 + 1] += f12.x[1]; g->d[ff->i1*3 + 2] += f12.x[2]; g->d[ff->i2*3 + 0] -= f12.x[0]; g->d[ff->i2*3 + 1] -= f12.x[1]; g->d[ff->i2*3 + 2] -= f12.x[2]; g->d[ff->i3*3 + 0] -= f34.x[0]; g->d[ff->i3*3 + 1] -= f34.x[1]; g->d[ff->i3*3 + 2] -= f34.x[2]; g->d[ff->i4*3 + 0] += f34.x[0]; g->d[ff->i4*3 + 1] += f34.x[1]; g->d[ff->i4*3 + 2] += f34.x[2]; } ee = 0.5*d*d*ff->fc; e->d[0] += ee; printf("Restraint %2d bd2 : value %11.4e energy %11.4e force %11.4e\n",ff->index,r12-r34,ee,-d*ff->fc); break; case BONDDIFF3: a1 = &f->atoms[ff->i1]; a2 = &f->atoms[ff->i2]; a3 = &f->atoms[ff->i3]; a4 = &f->atoms[ff->i4]; a5 = &f->atoms[ff->i5]; a6 = &f->atoms[ff->i6]; r12 = vdist(a1->pos,a2->pos); r34 = vdist(a3->pos,a4->pos); r56 = vdist(a5->pos,a6->pos); v12 = vnorm(vdiff(a2->pos,a1->pos)); v34 = vnorm(vdiff(a4->pos,a3->pos)); v56 = vnorm(vdiff(a6->pos,a5->pos)); d = r12 - r34 -r56 - ff->r0; f12 = vscal(v12,-d*ff->fc); f34 = vscal(v34,-d*ff->fc); f56 = vscal(v56,-d*ff->fc); if(g->d){ g->d[ff->i1*3 + 0] += f12.x[0]; g->d[ff->i1*3 + 1] += f12.x[1]; g->d[ff->i1*3 + 2] += f12.x[2]; g->d[ff->i2*3 + 0] -= f12.x[0]; g->d[ff->i2*3 + 1] -= f12.x[1]; g->d[ff->i2*3 + 2] -= f12.x[2]; g->d[ff->i3*3 + 0] -= f34.x[0]; g->d[ff->i3*3 + 1] -= f34.x[1]; g->d[ff->i3*3 + 2] -= f34.x[2]; g->d[ff->i4*3 + 0] += f34.x[0]; g->d[ff->i4*3 + 1] += f34.x[1]; g->d[ff->i4*3 + 2] += f34.x[2]; g->d[ff->i5*3 + 0] -= f56.x[0]; g->d[ff->i5*3 + 1] -= f56.x[1]; g->d[ff->i5*3 + 2] -= f56.x[2]; g->d[ff->i6*3 + 0] += f56.x[0]; g->d[ff->i6*3 + 1] += f56.x[1]; g->d[ff->i6*3 + 2] += f56.x[2]; } ee = 0.5*d*d*ff->fc; e->d[0] += ee; printf("Restraint %2d bd3 : value %11.4e energy %11.4e force %11.4e\n",ff->index,r12-r34-r56,ee,-d*ff->fc); break; case BONDDIFF4: a1 = &f->atoms[ff->i1]; a2 = &f->atoms[ff->i2]; a3 = &f->atoms[ff->i3]; a4 = &f->atoms[ff->i4]; a5 = &f->atoms[ff->i5]; a6 = &f->atoms[ff->i6]; a7 = &f->atoms[ff->i7]; a8 = &f->atoms[ff->i8]; r12 = vdist(a1->pos,a2->pos); r34 = vdist(a3->pos,a4->pos); r56 = vdist(a5->pos,a6->pos); r78 = vdist(a7->pos,a8->pos); v12 = vnorm(vdiff(a2->pos,a1->pos)); v34 = vnorm(vdiff(a4->pos,a3->pos)); v56 = vnorm(vdiff(a6->pos,a5->pos)); v78 = vnorm(vdiff(a8->pos,a7->pos)); d = r12 + r34 - r56 - r78 - ff->r0; f12 = vscal(v12,-d*ff->fc); f34 = vscal(v34,-d*ff->fc); f56 = vscal(v56,-d*ff->fc); f78 = vscal(v78,-d*ff->fc); if(g->d){ g->d[ff->i1*3 + 0] += f12.x[0]; g->d[ff->i1*3 + 1] += f12.x[1]; g->d[ff->i1*3 + 2] += f12.x[2]; g->d[ff->i2*3 + 0] -= f12.x[0]; g->d[ff->i2*3 + 1] -= f12.x[1]; g->d[ff->i2*3 + 2] -= f12.x[2]; g->d[ff->i3*3 + 0] += f34.x[0]; g->d[ff->i3*3 + 1] += f34.x[1]; g->d[ff->i3*3 + 2] += f34.x[2]; g->d[ff->i4*3 + 0] -= f34.x[0]; g->d[ff->i4*3 + 1] -= f34.x[1]; g->d[ff->i4*3 + 2] -= f34.x[2]; g->d[ff->i5*3 + 0] -= f56.x[0]; g->d[ff->i5*3 + 1] -= f56.x[1]; g->d[ff->i5*3 + 2] -= f56.x[2]; g->d[ff->i6*3 + 0] += f56.x[0]; g->d[ff->i6*3 + 1] += f56.x[1]; g->d[ff->i6*3 + 2] += f56.x[2]; g->d[ff->i7*3 + 0] -= f78.x[0]; g->d[ff->i7*3 + 1] -= f78.x[1]; g->d[ff->i7*3 + 2] -= f78.x[2]; g->d[ff->i8*3 + 0] += f78.x[0]; g->d[ff->i8*3 + 1] += f78.x[1]; g->d[ff->i8*3 + 2] += f78.x[2]; } ee = 0.5*d*d*ff->fc; e->d[0] += ee; printf("Restraint %2d bd4 : value %11.4e energy %11.4e force %11.4e\n",ff->index,r12+r34-r56-r78,ee,-d*ff->fc); break; case SBC_ATOM: ee= 0.; r32=0.; for(iat=0; iat < f->natom; iat++){ a1 = &f->atoms[ff->i1]; a2 = &f->atoms[iat]; /* Centre */ r12 = vdist(a1->pos,a2->pos); if(r12 > r32) r32=r12; d = r12 - ff->r0; if(d > 0.0) { v12 = vnorm(vdiff(a2->pos,a1->pos)); f12 = vscal(v12,d*ff->fc); if(g->d){ g->d[iat*3 + 0] += f12.x[0]; g->d[iat*3 + 1] += f12.x[1]; g->d[iat*3 + 2] += f12.x[2]; g->d[ff->i1 + 0] -= f12.x[0]; g->d[ff->i1 + 1] -= f12.x[1]; g->d[ff->i1 + 2] -= f12.x[2]; } ee += 0.5*d*d*ff->fc; } } e->d[0] += ee; printf("Restraint %2d SBC-A: value %11.4e energy %11.4e\n",ff->index,r32,ee); break; case SBC: ee= 0.; r32=0.; pressure=0.; pi = 4. * atan(1.); for(iat=0; iat < f->natom; iat++){ b1 = f->atoms[iat]; b2.pos.x[0]= ff->x; b2.pos.x[1]= ff->y; b2.pos.x[2]= ff->z; r12 = vdist(b1.pos,b2.pos); if(r12 > r32) r32=r12; d = r12 - ff->r0; /* ignore the center atom */ if(d > 0.0) { v12 = vnorm(vdiff(b2.pos,b1.pos)); f12 = vscal(v12,d*ff->fc); if(g->d){ g->d[iat*3 + 0] -= f12.x[0]; g->d[iat*3 + 1] -= f12.x[1]; g->d[iat*3 + 2] -= f12.x[2]; } pressure += d * ff->fc / (4.*pi*r12*r12); ee += 0.5*d*d*ff->fc; } } e->d[0] += ee; printf("Restraint %2d SBC : value %11.4e energy %11.4e\n",ff->index,r32,ee); /* printf("SBC pressure: %11.4e a.u. = %11.4e bar\n",pressure, pressure * 2.9421011987E+8 ); */ break; case SBC_IND: ee= 0.; r32=0.; pressure=0.; pi = 4. * atan(1.); b1 = f->atoms[ff->i1]; b2.pos.x[0]= ff->x; b2.pos.x[1]= ff->y; b2.pos.x[2]= ff->z; r12 = vdist(b1.pos,b2.pos); if(r12 > r32) r32=r12; d = r12 - ff->r0; /* ignore the center atom */ if(d > 0.0) { v12 = vnorm(vdiff(b2.pos,b1.pos)); f12 = vscal(v12,d*ff->fc); if(g->d){ g->d[ff->i1*3 + 0] -= f12.x[0]; g->d[ff->i1*3 + 1] -= f12.x[1]; g->d[ff->i1*3 + 2] -= f12.x[2]; } pressure += d * ff->fc / (4.*pi*r12*r12); ee = 0.5*d*d*ff->fc; } e->d[0] += ee; if(debug)printf("Restraint %2d SBC_individual : value %11.4e energy %11.4e\n",ff->index,r32,ee); break; case CART: /* handle cart and cart_atom on an equal basis */ case CART_ATOM: b1 = f->atoms[ff->i1]; b2.pos.x[0]= ff->x; b2.pos.x[1]= ff->y; b2.pos.x[2]= ff->z; r12 = vdist(b1.pos,b2.pos); d = r12; if(r12 > 0. ) { v12 = vnorm(vdiff(b2.pos,b1.pos)); f12 = vscal(v12,-d*ff->fc); if(g->d){ g->d[ff->i1*3 + 0] += f12.x[0]; g->d[ff->i1*3 + 1] += f12.x[1]; g->d[ff->i1*3 + 2] += f12.x[2]; } } ee = 0.5*d*d*ff->fc; e->d[0] += ee; if(debug)printf("r,energy in calc_ff %lf %10.4e\n",r12,ee); printf("Restraint %2d cart : value %11.4e energy %11.4e force %11.4e\n",ff->index,r12,ee,-d*ff->fc); break; default: printf("energy term not implemented\n"); exit(1); } fflush(stdout); return 0; } /* new_ff structure, to add an ff maintaining order */ Ff new_rst_ff(struct rst_struct *rst) { Ff ff; ff = New(ff_s); if(rst->current_ff){ rst->current_ff->prev = ff; }else{ rst->last_ff = ff; } rst->current_ff = ff; ff->prev = NULL; return ff; } /* Here starts the Block for the dispersion correction following Grimme */ /* S.Grimme, J.Comput.Chem. 27, (2006), 1787-1799. */ int DispersionCmd(ClientData clientData, Tcl_Interp *interp, int argc, const char *argv[]) { /* Local variables */ Frag coords; Matrix energy, gradient; ObjList lcoords, lene, lgrad; const char *s; int nat,ipar,iat,jat,iz,jz; /* string-buffers for use in Tcl_Eval , Tcl_SetVar (to avoid conflict with strbuff altered inside the called C-routines */ char tclbuff[STRBUFFLEN]; char tclbuff2[STRBUFFLEN]; char buff2[2]; #define MAXPAR 100 double rnod[MAXPAR],c6[MAXPAR]; double rsum,dist,c6ij,dexp; double edisp,fdamp,s6,grad,edamp,c_over_r6; struct atom_struct *a1, *a2; struct vector_struct v12, f12; int verbose; /* how much output to generate: 0 or 1 */ int first; /* if 1 do not calculate, but test if parameters are set, write citation info */ /* scaling factor for the funcational */ GetVar(s,"s6"); if(Tcl_GetDouble(interp,s,&s6)){ printf("bad s6 in dispersion: %s\n",s); return TCL_ERROR; } /* verbosity */ GetVar(s,"verbose"); if(Tcl_GetInt(interp,s,&verbose)){ printf("bad verbose in dispersion: %s\n",s); return TCL_ERROR; } GetVar(s,"first"); if(Tcl_GetInt(interp,s,&first)){ printf("bad first in dispersion: %s\n",s); return TCL_ERROR; } /* Access initial structure */ GetVar(s,"coords"); lcoords = get_objlist(s,"fragment",CHEMSH_OBJ_EXISTS,0); if(!lcoords){ printf("bad coords tag in dispersion_correction: %s\n",s); return TCL_ERROR; } coords = (Frag) lcoords->object->data; nat = FRAG_get_number_of_atoms(coords); /* Gradient matrix to accumulate into (3n*1) */ GetVar(s,"gradient"); lgrad = get_objlist(s,"matrix",CHEMSH_OBJ_EXISTS,0); if(!lgrad){ /* printf("bad matrix tag for gradient in dispersion_correction: %s \n",s); return TCL_ERROR; */ gradient = NULL; } else { gradient = (Matrix) lgrad->object->data; } /* energy matrix (1*1) */ GetVar(s,"energy"); lene = get_objlist(s,"matrix",CHEMSH_OBJ_EXISTS,0); if(!lene){ printf("bad matrix tag for energy in dispersion_correction: %s\n", s); return TCL_ERROR; } energy = (Matrix) lene->object->data; /* set parameters according to Table 1 in S.Grimme, J.Comput.Chem. 27, (2006), 1787-1799. rnod is R0, the vdW radius c6 is the C6 coefficient per atom rsum is the summ of the vdW radii of the two atoms dist is the distance between the atoms c6ij is the geometric mean of the c6 parameters of the two atoms dexp is the exponential factor d in Eq. 12 in Grimme's paper */ /* set default values for atom types that are not parametrised */ for(ipar=1;iparatoms[iat]; iz = a1->znum; if(iz==0){ printf("Error: could not determine Z of atom %d\n", iat+1); return TCL_ERROR; } if(iz>100){ printf("Error: dispersion correction for atom %d (Element %d) not possible!\n", iat+1,iz); return TCL_ERROR; } if(iz>54 && first){ printf("Warning: dispersion correction for atom %d (Element %d) not parametrised!\n", iat+1,iz); } for(jat=iat+1;jatatoms[jat]; jz= a2->znum; if(jz==0){ printf("Error: j could not determine Z of atom %d\n", jat+1); return TCL_ERROR; } if(jz>100){ printf("Error: dispersion correction for atom %d (Element %d) not parametrised!\n", jat+1,jz); return TCL_ERROR; } if(jz>54 && first){ printf("Warning: dispersion correction for atom %d (Element %d) not implemented!\n", jat+1,jz); } /* Implementation of Equations 11, 12, and 13 of the Grimme paper */ rsum = (rnod[iz] + rnod[jz] ) / 0.52917706 ; /* convert ang to au */ c6ij=sqrt( c6[iz] * c6[jz] ) * 17.345304 ; /* convert J*nm^6/mol to au */ a1 = &coords->atoms[iat]; a2 = &coords->atoms[jat]; dist = vdist(a1->pos,a2->pos); v12 = vnorm(vdiff(a2->pos,a1->pos)); edamp = exp( -dexp * (dist/rsum - 1.) ) ; c_over_r6 = c6ij / ( dist *dist *dist *dist *dist *dist); fdamp=1./( 1. + edamp ); edisp += -s6 * c_over_r6 * fdamp; if(gradient) { if(gradient->d && !first){ grad = s6 * ( -6. * c_over_r6 * fdamp / dist + c_over_r6 * dexp * edamp / rsum / ( 1. + edamp ) / ( 1. + edamp ) ); f12 = vscal(v12,grad); gradient->d[iat*3 + 0] += f12.x[0]; gradient->d[iat*3 + 1] += f12.x[1]; gradient->d[iat*3 + 2] += f12.x[2]; gradient->d[jat*3 + 0] -= f12.x[0]; gradient->d[jat*3 + 1] -= f12.x[1]; gradient->d[jat*3 + 2] -= f12.x[2]; } } } } if( ! first ) { energy->d[0] += edisp ; if(verbose >= 1) printf("Dispersion correction included in the energy: %10.5f Hartree\n",edisp); } rel_objlist(lcoords); rel_objlist(lene); if(gradient) rel_objlist(lgrad); return TCL_OK; }