TclCommands.C File Reference

#include <stdlib.h>
#include <malloc.h>
#include <errno.h>
#include "TclCommands.h"
#include "Vector.h"
#include "NamdTypes.h"
#include <tcl.h>
#include "Matrix4.C"
#include "TclVec.C"

Go to the source code of this file.

Defines

#define USE_COMPAT_CONST

Functions

static int get_3D_vector (Tcl_Interp *interp, Tcl_Obj *const list, Vector &result)
static Tcl_Obj * obj_3D_vector (const Vector &v)
int proc_getbond (ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
int proc_getangle (ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
int proc_getdihedral (ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
int proc_anglegrad (ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
int proc_dihedralgrad (ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
int tcl_vector_math_init (Tcl_Interp *interp)


Define Documentation

#define USE_COMPAT_CONST

Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved.

Definition at line 18 of file TclCommands.C.


Function Documentation

static int get_3D_vector ( Tcl_Interp *  interp,
Tcl_Obj *const   list,
Vector result 
) [static]

Definition at line 26 of file TclCommands.C.

References Vector::x, Vector::y, and Vector::z.

Referenced by proc_anglegrad(), proc_dihedralgrad(), proc_getangle(), proc_getbond(), and proc_getdihedral().

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 }

static Tcl_Obj* obj_3D_vector ( const Vector v  )  [static]

Definition at line 45 of file TclCommands.C.

References Vector::x, Vector::y, and Vector::z.

Referenced by proc_anglegrad(), and proc_dihedralgrad().

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 }

int proc_anglegrad ( ClientData  ,
Tcl_Interp *  interp,
int  argc,
Tcl_Obj *const   argv[] 
)

Definition at line 131 of file TclCommands.C.

References get_3D_vector(), Vector::length(), and obj_3D_vector().

Referenced by tcl_vector_math_init().

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 }

int proc_dihedralgrad ( ClientData  ,
Tcl_Interp *  interp,
int  argc,
Tcl_Obj *const   argv[] 
)

Definition at line 181 of file TclCommands.C.

References A, B, cross(), get_3D_vector(), Vector::length(), obj_3D_vector(), Vector::x, Vector::y, and Vector::z.

Referenced by tcl_vector_math_init().

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 }

int proc_getangle ( ClientData  ,
Tcl_Interp *  interp,
int  argc,
Tcl_Obj *const   argv[] 
)

Definition at line 76 of file TclCommands.C.

References get_3D_vector(), Vector::length(), and PI.

Referenced by tcl_vector_math_init().

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 }

int proc_getbond ( ClientData  ,
Tcl_Interp *  interp,
int  argc,
Tcl_Obj *const   argv[] 
)

Definition at line 59 of file TclCommands.C.

References get_3D_vector().

Referenced by tcl_vector_math_init().

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 }

int proc_getdihedral ( ClientData  ,
Tcl_Interp *  interp,
int  argc,
Tcl_Obj *const   argv[] 
)

Definition at line 98 of file TclCommands.C.

References A, B, cross(), get_3D_vector(), Vector::length(), and PI.

Referenced by tcl_vector_math_init().

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 }

int tcl_vector_math_init ( Tcl_Interp *  interp  ) 

Definition at line 300 of file TclCommands.C.

References proc_anglegrad(), proc_dihedralgrad(), proc_getangle(), proc_getbond(), proc_getdihedral(), and Vec_Init().

Referenced by ComputeTclBC::ComputeTclBC(), ScriptTcl::ScriptTcl(), and ScriptTcl::tclsh().

00300                                              {
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 }


Generated on Sat Sep 23 01:17:17 2017 for NAMD by  doxygen 1.4.7