TclCommands.C

Go to the documentation of this file.
00001 
00007 #include <stdlib.h> 
00008 #ifndef _NO_MALLOC_H
00009 #include <malloc.h>
00010 #endif
00011 #include <errno.h>
00012 #include "TclCommands.h"
00013 #include "Vector.h"
00014 #include "NamdTypes.h"
00015 
00016 #ifdef NAMD_TCL
00017 
00018 #define USE_COMPAT_CONST
00019 #include <tcl.h>
00020 
00021 #include "Matrix4.C"
00022 #include "TclVec.C"
00023 
00024 
00025 // Get a 3-D vector from a TCL list
00026 static int get_3D_vector(Tcl_Interp *interp, Tcl_Obj *const list, Vector &result)
00027 {
00028   int num;
00029   Tcl_Obj **data;
00030   
00031   if (Tcl_ListObjGetElements(interp, list, &num, &data) != TCL_OK)
00032     return 0;
00033   if (Tcl_GetDoubleFromObj(interp,data[0],&(result.x)) != TCL_OK)
00034     return 0;
00035   if (Tcl_GetDoubleFromObj(interp,data[1],&(result.y)) != TCL_OK)
00036     return 0;
00037   if (Tcl_GetDoubleFromObj(interp,data[2],&(result.z)) != TCL_OK)
00038     return 0;
00039 
00040   return 1;
00041 }
00042 
00043 
00044 // Append a 3-D vector to the result string
00045 static Tcl_Obj* obj_3D_vector(const Vector &v)
00046 {
00047   Tcl_Obj* doublev[3];
00048   
00049   doublev[0] = Tcl_NewDoubleObj(v.x);
00050   doublev[1] = Tcl_NewDoubleObj(v.y);
00051   doublev[2] = Tcl_NewDoubleObj(v.z);
00052   
00053   return Tcl_NewListObj(3,doublev);
00054 }
00055 
00056 
00057 // Function: getbond coor1 coor2
00058 //  Returns: the length of the bond formed by the two atoms (i.e., the distance between them)
00059 int proc_getbond(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
00060 {
00061   Vector r1, r2;
00062   
00063   if (argc != 3)
00064     return TCL_ERROR;
00065   if (!get_3D_vector(interp,argv[1],r1))
00066     return TCL_ERROR;
00067   if (!get_3D_vector(interp,argv[2],r2))
00068     return TCL_ERROR;
00069   Tcl_SetObjResult(interp, Tcl_NewDoubleObj((r2-r1).length()));
00070   return TCL_OK;
00071 }
00072 
00073 
00074 // Function: getangle coor1 coor2 coor3
00075 //  Returns: the angle formed by the three atoms
00076 int proc_getangle(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
00077 {
00078   Vector r1, r2, r3, r12, r32;
00079   
00080   if (argc != 4)
00081     return TCL_ERROR;
00082   if (!get_3D_vector(interp,argv[1],r1))
00083     return TCL_ERROR;
00084   if (!get_3D_vector(interp,argv[2],r2))
00085     return TCL_ERROR;
00086   if (!get_3D_vector(interp,argv[3],r3))
00087     return TCL_ERROR;
00088   r12 = r1 - r2;
00089   r32 = r3 - r2;
00090   Tcl_SetObjResult(interp, Tcl_NewDoubleObj(
00091                 acos((r12*r32)/(r12.length()*r32.length()))*180/PI  ));
00092   return TCL_OK;
00093 }
00094 
00095 
00096 // Function: getdihedral coor1 coor2 coor3 coor4
00097 //  Returns: the dihedral formed by the four atoms
00098 int proc_getdihedral(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
00099 {
00100   BigReal rA, rB, rC;
00101   Vector r1, r2, r3, r4, r12, r23, r34, A, B, C;
00102   
00103   if (argc != 5)
00104     return TCL_ERROR;
00105   if (!get_3D_vector(interp,argv[1],r1))
00106     return TCL_ERROR;
00107   if (!get_3D_vector(interp,argv[2],r2))
00108     return TCL_ERROR;
00109   if (!get_3D_vector(interp,argv[3],r3))
00110     return TCL_ERROR;
00111   if (!get_3D_vector(interp,argv[4],r4))
00112     return TCL_ERROR;
00113   r12 = r1 - r2;
00114   r23 = r2 - r3;
00115   r34 = r3 - r4;
00116   A = cross(r12,r23);
00117   B = cross(r23,r34);
00118   C = cross(r23,A);
00119   rA = A.length();
00120   rB = B.length();
00121   rC = C.length();
00122   Tcl_SetObjResult(interp, Tcl_NewDoubleObj(
00123                 -atan2((C*B)/(rC*rB),(A*B)/(rA*rB))*180/PI  ));
00124   return TCL_OK;
00125 }
00126 
00127 
00128 // Function: anglegrad coor1 coor2 coor3
00129 //  Returns: a list of gradients for each atom
00130 // The code was basically copied from ComputeAngles.C
00131 int proc_anglegrad(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
00132 {
00133   Vector r1, r2, r3, r12, r32;
00134   
00135   if (argc != 4)
00136     return TCL_ERROR;
00137   if (!get_3D_vector(interp,argv[1],r1))
00138     return TCL_ERROR;
00139   if (!get_3D_vector(interp,argv[2],r2))
00140     return TCL_ERROR;
00141   if (!get_3D_vector(interp,argv[3],r3))
00142     return TCL_ERROR;
00143   
00144   r12 = r1 - r2;
00145   BigReal d12 = r12.length();
00146   r32 = r3 - r2;
00147   BigReal d32 = r32.length();
00148   
00149   BigReal cos_theta = (r12*r32)/(d12*d32);
00150   
00151   //  Normalize vector r12 and r32
00152   BigReal d12inv = 1. / d12;
00153   BigReal d32inv = 1. / d32;
00154   
00155   BigReal sin_theta = sqrt(1.0 - cos_theta*cos_theta);
00156   BigReal diff = 1/sin_theta;
00157   BigReal c1 = diff * d12inv;
00158   BigReal c2 = diff * d32inv;
00159   
00160   //  Calculate the actual forces
00161   Force force1 = c1*(r12*(d12inv*cos_theta) - r32*d32inv);
00162   Force force2 = force1;
00163   Force force3 = c2*(r32*(d32inv*cos_theta) - r12*d12inv);
00164   force2 += force3;  force2 *= -1;
00165   
00166   Tcl_Obj* forcev[3];
00167 
00168   forcev[0] = obj_3D_vector(force1);
00169   forcev[1] = obj_3D_vector(force2);
00170   forcev[2] = obj_3D_vector(force3);
00171 
00172   Tcl_SetObjResult(interp, Tcl_NewListObj(3, forcev));
00173   
00174   return TCL_OK;
00175 }
00176 
00177 
00178 // Function: dihedralgrad coor1 coor2 coor3 coor4
00179 //  Returns: a list of gradients for each atom
00180 // The code was basically copied from ComputeDihedrals.C
00181 int proc_dihedralgrad(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
00182 {
00183   BigReal K1;
00184   Vector r1, r2, r3, r4, r12, r23, r34;
00185   
00186   if (argc != 5)
00187     return TCL_ERROR;
00188   if (!get_3D_vector(interp,argv[1],r1))
00189     return TCL_ERROR;
00190   if (!get_3D_vector(interp,argv[2],r2))
00191     return TCL_ERROR;
00192   if (!get_3D_vector(interp,argv[3],r3))
00193     return TCL_ERROR;
00194   if (!get_3D_vector(interp,argv[4],r4))
00195     return TCL_ERROR;
00196   
00197   r12 = r1 - r2;
00198   r23 = r2 - r3;
00199   r34 = r3 - r4;
00200   
00201   //  Calculate the cross products and distances
00202   Vector A = cross(r12,r23);
00203   BigReal rA = A.length();
00204   Vector B = cross(r23,r34);
00205   BigReal rB = B.length();
00206   Vector C = cross(r23,A);
00207   BigReal rC = C.length();
00208   
00209   //  Calculate the sin and cos
00210   BigReal cos_phi = (A*B)/(rA*rB);
00211   BigReal sin_phi = (C*B)/(rC*rB);
00212   
00213   Force f1,f2,f3;
00214   
00215   //  Normalize B
00216   rB = 1.0/rB;
00217   B *= rB;
00218   
00219   //  We first need to figure out whether the
00220   //  sin or cos form will be more stable.  For this,
00221   //  just look at the value of phi
00222   if (fabs(sin_phi) > 0.1)
00223   {
00224     //  use the sin version to avoid 1/cos terms
00225     
00226     //  Normalize A
00227     rA = 1.0/rA;
00228     A *= rA;
00229     Vector dcosdA = rA*(cos_phi*A-B);
00230     Vector dcosdB = rB*(cos_phi*B-A);
00231     
00232     K1 = -1/sin_phi;
00233     
00234     f1.x = K1*(r23.y*dcosdA.z - r23.z*dcosdA.y);
00235     f1.y = K1*(r23.z*dcosdA.x - r23.x*dcosdA.z);
00236     f1.z = K1*(r23.x*dcosdA.y - r23.y*dcosdA.x);
00237     
00238     f3.x = K1*(r23.z*dcosdB.y - r23.y*dcosdB.z);
00239     f3.y = K1*(r23.x*dcosdB.z - r23.z*dcosdB.x);
00240     f3.z = K1*(r23.y*dcosdB.x - r23.x*dcosdB.y);
00241     
00242     f2.x = K1*(r12.z*dcosdA.y - r12.y*dcosdA.z
00243              + r34.y*dcosdB.z - r34.z*dcosdB.y);
00244     f2.y = K1*(r12.x*dcosdA.z - r12.z*dcosdA.x
00245              + r34.z*dcosdB.x - r34.x*dcosdB.z);
00246     f2.z = K1*(r12.y*dcosdA.x - r12.x*dcosdA.y
00247              + r34.x*dcosdB.y - r34.y*dcosdB.x);
00248   }
00249   else
00250   {
00251     //  This angle is closer to 0 or 180 than it is to
00252     //  90, so use the cos version to avoid 1/sin terms
00253     
00254     //  Normalize C
00255     rC = 1.0/rC;
00256     C *= rC;
00257     Vector dsindC = rC*(sin_phi*C-B);
00258     Vector dsindB = rB*(sin_phi*B-C);
00259     
00260     K1 = 1/cos_phi;
00261     
00262     f1.x = K1*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
00263               - r23.x*r23.y*dsindC.y
00264               - r23.x*r23.z*dsindC.z);
00265     f1.y = K1*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
00266               - r23.y*r23.z*dsindC.z
00267               - r23.y*r23.x*dsindC.x);
00268     f1.z = K1*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
00269               - r23.z*r23.x*dsindC.x
00270               - r23.z*r23.y*dsindC.y);
00271     
00272     f3 = cross(K1,dsindB,r23);
00273     
00274     f2.x = K1*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
00275            +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
00276            +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
00277            +dsindB.z*r34.y - dsindB.y*r34.z);
00278     f2.y = K1*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
00279            +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
00280            +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
00281            +dsindB.x*r34.z - dsindB.z*r34.x);
00282     f2.z = K1*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
00283            +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
00284            +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
00285            +dsindB.y*r34.x - dsindB.x*r34.y);
00286   }
00287   
00288   Tcl_Obj* forcev[4];
00289 
00290   forcev[0] = obj_3D_vector(f1);
00291   forcev[1] = obj_3D_vector(f2-f1);
00292   forcev[2] = obj_3D_vector(f3-f2);
00293   forcev[3] = obj_3D_vector(-f3);
00294 
00295   Tcl_SetObjResult(interp, Tcl_NewListObj(4, forcev));
00296   
00297   return TCL_OK;
00298 }
00299 
00300 int tcl_vector_math_init(Tcl_Interp *interp) {
00301 
00302   // first import from TclVec.C stolen from VMD
00303   Vec_Init(interp);
00304 
00305   Tcl_CreateObjCommand(interp, "getbond", proc_getbond,
00306     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00307   Tcl_CreateObjCommand(interp, "getangle", proc_getangle,
00308     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00309   Tcl_CreateObjCommand(interp, "getdihedral", proc_getdihedral,
00310     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00311   Tcl_CreateObjCommand(interp, "anglegrad", proc_anglegrad,
00312     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00313   Tcl_CreateObjCommand(interp, "dihedralgrad", proc_dihedralgrad,
00314     (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
00315 
00316   return TCL_OK;
00317 }
00318 
00319 
00320 #endif

Generated on Thu Nov 23 01:17:14 2017 for NAMD by  doxygen 1.4.7