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.

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)

Function Documentation

get_3D_vector()

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

static

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

Definition at line 25 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().

{
  int num;
  Tcl_Obj **data;
  
  if (Tcl_ListObjGetElements(interp, list, &num, &data) != TCL_OK)
    return 0;
  if (Tcl_GetDoubleFromObj(interp,data[0],&(result.x)) != TCL_OK)
    return 0;
  if (Tcl_GetDoubleFromObj(interp,data[1],&(result.y)) != TCL_OK)
    return 0;
  if (Tcl_GetDoubleFromObj(interp,data[2],&(result.z)) != TCL_OK)
    return 0;

  return 1;
}

obj_3D_vector()

static Tcl_Obj* obj_3D_vector ( const Vector &  v )

static

Definition at line 44 of file TclCommands.C.

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

Referenced by proc_anglegrad(), and proc_dihedralgrad().

{
  Tcl_Obj* doublev[3];
  
  doublev[0] = Tcl_NewDoubleObj(v.x);
  doublev[1] = Tcl_NewDoubleObj(v.y);
  doublev[2] = Tcl_NewDoubleObj(v.z);
  
  return Tcl_NewListObj(3,doublev);
}

proc_anglegrad()

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

Definition at line 130 of file TclCommands.C.

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

Referenced by tcl_vector_math_init().

{
  Vector r1, r2, r3, r12, r32;
  
  if (argc != 4)
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[1],r1))
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[2],r2))
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[3],r3))
    return TCL_ERROR;
  
  r12 = r1 - r2;
  BigReal d12 = r12.length();
  r32 = r3 - r2;
  BigReal d32 = r32.length();
  
  BigReal cos_theta = (r12*r32)/(d12*d32);
  
  //  Normalize vector r12 and r32
  BigReal d12inv = 1. / d12;
  BigReal d32inv = 1. / d32;
  
  BigReal sin_theta = sqrt(1.0 - cos_theta*cos_theta);
  BigReal diff = 1/sin_theta;
  BigReal c1 = diff * d12inv;
  BigReal c2 = diff * d32inv;
  
  //  Calculate the actual forces
  Force force1 = c1*(r12*(d12inv*cos_theta) - r32*d32inv);
  Force force2 = force1;
  Force force3 = c2*(r32*(d32inv*cos_theta) - r12*d12inv);
  force2 += force3;  force2 *= -1;
  
  Tcl_Obj* forcev[3];

  forcev[0] = obj_3D_vector(force1);
  forcev[1] = obj_3D_vector(force2);
  forcev[2] = obj_3D_vector(force3);

  Tcl_SetObjResult(interp, Tcl_NewListObj(3, forcev));
  
  return TCL_OK;
}

proc_dihedralgrad()

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

Definition at line 180 of file TclCommands.C.

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

Referenced by tcl_vector_math_init().

{
  BigReal K1;
  Vector r1, r2, r3, r4, r12, r23, r34;
  
  if (argc != 5)
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[1],r1))
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[2],r2))
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[3],r3))
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[4],r4))
    return TCL_ERROR;
  
  r12 = r1 - r2;
  r23 = r2 - r3;
  r34 = r3 - r4;
  
  //  Calculate the cross products and distances
  Vector A = cross(r12,r23);
  BigReal rA = A.length();
  Vector B = cross(r23,r34);
  BigReal rB = B.length();
  Vector C = cross(r23,A);
  BigReal rC = C.length();
  
  //  Calculate the sin and cos
  BigReal cos_phi = (A*B)/(rA*rB);
  BigReal sin_phi = (C*B)/(rC*rB);
  
  Force f1,f2,f3;
  
  //  Normalize B
  rB = 1.0/rB;
  B *= rB;
  
  //  We first need to figure out whether the
  //  sin or cos form will be more stable.  For this,
  //  just look at the value of phi
  if (fabs(sin_phi) > 0.1)
  {
    //  use the sin version to avoid 1/cos terms
    
    //  Normalize A
    rA = 1.0/rA;
    A *= rA;
    Vector dcosdA = rA*(cos_phi*A-B);
    Vector dcosdB = rB*(cos_phi*B-A);
    
    K1 = -1/sin_phi;
    
    f1.x = K1*(r23.y*dcosdA.z - r23.z*dcosdA.y);
    f1.y = K1*(r23.z*dcosdA.x - r23.x*dcosdA.z);
    f1.z = K1*(r23.x*dcosdA.y - r23.y*dcosdA.x);
    
    f3.x = K1*(r23.z*dcosdB.y - r23.y*dcosdB.z);
    f3.y = K1*(r23.x*dcosdB.z - r23.z*dcosdB.x);
    f3.z = K1*(r23.y*dcosdB.x - r23.x*dcosdB.y);
    
    f2.x = K1*(r12.z*dcosdA.y - r12.y*dcosdA.z
             + r34.y*dcosdB.z - r34.z*dcosdB.y);
    f2.y = K1*(r12.x*dcosdA.z - r12.z*dcosdA.x
             + r34.z*dcosdB.x - r34.x*dcosdB.z);
    f2.z = K1*(r12.y*dcosdA.x - r12.x*dcosdA.y
             + r34.x*dcosdB.y - r34.y*dcosdB.x);
  }
  else
  {
    //  This angle is closer to 0 or 180 than it is to
    //  90, so use the cos version to avoid 1/sin terms
    
    //  Normalize C
    rC = 1.0/rC;
    C *= rC;
    Vector dsindC = rC*(sin_phi*C-B);
    Vector dsindB = rB*(sin_phi*B-C);
    
    K1 = 1/cos_phi;
    
    f1.x = K1*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
              - r23.x*r23.y*dsindC.y
              - r23.x*r23.z*dsindC.z);
    f1.y = K1*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
              - r23.y*r23.z*dsindC.z
              - r23.y*r23.x*dsindC.x);
    f1.z = K1*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
              - r23.z*r23.x*dsindC.x
              - r23.z*r23.y*dsindC.y);
    
    f3 = cross(K1,dsindB,r23);
    
    f2.x = K1*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
           +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
           +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
           +dsindB.z*r34.y - dsindB.y*r34.z);
    f2.y = K1*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
           +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
           +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
           +dsindB.x*r34.z - dsindB.z*r34.x);
    f2.z = K1*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
           +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
           +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
           +dsindB.y*r34.x - dsindB.x*r34.y);
  }
  
  Tcl_Obj* forcev[4];

  forcev[0] = obj_3D_vector(f1);
  forcev[1] = obj_3D_vector(f2-f1);
  forcev[2] = obj_3D_vector(f3-f2);
  forcev[3] = obj_3D_vector(-f3);

  Tcl_SetObjResult(interp, Tcl_NewListObj(4, forcev));
  
  return TCL_OK;
}

proc_getangle()

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

Definition at line 75 of file TclCommands.C.

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

Referenced by tcl_vector_math_init().

{
  Vector r1, r2, r3, r12, r32;
  
  if (argc != 4)
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[1],r1))
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[2],r2))
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[3],r3))
    return TCL_ERROR;
  r12 = r1 - r2;
  r32 = r3 - r2;
  Tcl_SetObjResult(interp, Tcl_NewDoubleObj(
                acos((r12*r32)/(r12.length()*r32.length()))*180/PI  ));
  return TCL_OK;
}

proc_getbond()

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

Definition at line 58 of file TclCommands.C.

References get_3D_vector().

Referenced by tcl_vector_math_init().

{
  Vector r1, r2;
  
  if (argc != 3)
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[1],r1))
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[2],r2))
    return TCL_ERROR;
  Tcl_SetObjResult(interp, Tcl_NewDoubleObj((r2-r1).length()));
  return TCL_OK;
}

proc_getdihedral()

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

Definition at line 97 of file TclCommands.C.

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

Referenced by tcl_vector_math_init().

{
  BigReal rA, rB, rC;
  Vector r1, r2, r3, r4, r12, r23, r34, A, B, C;
  
  if (argc != 5)
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[1],r1))
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[2],r2))
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[3],r3))
    return TCL_ERROR;
  if (!get_3D_vector(interp,argv[4],r4))
    return TCL_ERROR;
  r12 = r1 - r2;
  r23 = r2 - r3;
  r34 = r3 - r4;
  A = cross(r12,r23);
  B = cross(r23,r34);
  C = cross(r23,A);
  rA = A.length();
  rB = B.length();
  rC = C.length();
  Tcl_SetObjResult(interp, Tcl_NewDoubleObj(
                -atan2((C*B)/(rC*rB),(A*B)/(rA*rB))*180/PI  ));
  return TCL_OK;
}

tcl_vector_math_init()

int tcl_vector_math_init

(

Tcl_Interp *

)

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

Definition at line 299 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().

                                             {

  // first import from TclVec.C stolen from VMD
  Vec_Init(interp);

  Tcl_CreateObjCommand(interp, "getbond", proc_getbond,
    (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
  Tcl_CreateObjCommand(interp, "getangle", proc_getangle,
    (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
  Tcl_CreateObjCommand(interp, "getdihedral", proc_getdihedral,
    (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
  Tcl_CreateObjCommand(interp, "anglegrad", proc_anglegrad,
    (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
  Tcl_CreateObjCommand(interp, "dihedralgrad", proc_dihedralgrad,
    (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);

  return TCL_OK;
}