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

Generated on Mon Aug 20 01:17:16 2018 for NAMD by  doxygen 1.4.7