NAMD
Functions
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

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().

26 {
27  int num;
28  Tcl_Obj **data;
29 
30  if (Tcl_ListObjGetElements(interp, list, &num, &data) != TCL_OK)
31  return 0;
32  if (Tcl_GetDoubleFromObj(interp,data[0],&(result.x)) != TCL_OK)
33  return 0;
34  if (Tcl_GetDoubleFromObj(interp,data[1],&(result.y)) != TCL_OK)
35  return 0;
36  if (Tcl_GetDoubleFromObj(interp,data[2],&(result.z)) != TCL_OK)
37  return 0;
38 
39  return 1;
40 }
BigReal z
Definition: Vector.h:66
BigReal x
Definition: Vector.h:66
BigReal y
Definition: Vector.h:66
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().

45 {
46  Tcl_Obj* doublev[3];
47 
48  doublev[0] = Tcl_NewDoubleObj(v.x);
49  doublev[1] = Tcl_NewDoubleObj(v.y);
50  doublev[2] = Tcl_NewDoubleObj(v.z);
51 
52  return Tcl_NewListObj(3,doublev);
53 }
BigReal z
Definition: Vector.h:66
BigReal x
Definition: Vector.h:66
BigReal y
Definition: Vector.h:66
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().

131 {
132  Vector r1, r2, r3, r12, r32;
133 
134  if (argc != 4)
135  return TCL_ERROR;
136  if (!get_3D_vector(interp,argv[1],r1))
137  return TCL_ERROR;
138  if (!get_3D_vector(interp,argv[2],r2))
139  return TCL_ERROR;
140  if (!get_3D_vector(interp,argv[3],r3))
141  return TCL_ERROR;
142 
143  r12 = r1 - r2;
144  BigReal d12 = r12.length();
145  r32 = r3 - r2;
146  BigReal d32 = r32.length();
147 
148  BigReal cos_theta = (r12*r32)/(d12*d32);
149 
150  // Normalize vector r12 and r32
151  BigReal d12inv = 1. / d12;
152  BigReal d32inv = 1. / d32;
153 
154  BigReal sin_theta = sqrt(1.0 - cos_theta*cos_theta);
155  BigReal diff = 1/sin_theta;
156  BigReal c1 = diff * d12inv;
157  BigReal c2 = diff * d32inv;
158 
159  // Calculate the actual forces
160  Force force1 = c1*(r12*(d12inv*cos_theta) - r32*d32inv);
161  Force force2 = force1;
162  Force force3 = c2*(r32*(d32inv*cos_theta) - r12*d12inv);
163  force2 += force3; force2 *= -1;
164 
165  Tcl_Obj* forcev[3];
166 
167  forcev[0] = obj_3D_vector(force1);
168  forcev[1] = obj_3D_vector(force2);
169  forcev[2] = obj_3D_vector(force3);
170 
171  Tcl_SetObjResult(interp, Tcl_NewListObj(3, forcev));
172 
173  return TCL_OK;
174 }
static int get_3D_vector(Tcl_Interp *interp, Tcl_Obj *const list, Vector &result)
Definition: TclCommands.C:25
Definition: Vector.h:64
BigReal length(void) const
Definition: Vector.h:169
static Tcl_Obj * obj_3D_vector(const Vector &v)
Definition: TclCommands.C:44
double BigReal
Definition: common.h:114
int proc_dihedralgrad ( ClientData  ,
Tcl_Interp *  interp,
int  argc,
Tcl_Obj *const  argv[] 
)

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

181 {
182  BigReal K1;
183  Vector r1, r2, r3, r4, r12, r23, r34;
184 
185  if (argc != 5)
186  return TCL_ERROR;
187  if (!get_3D_vector(interp,argv[1],r1))
188  return TCL_ERROR;
189  if (!get_3D_vector(interp,argv[2],r2))
190  return TCL_ERROR;
191  if (!get_3D_vector(interp,argv[3],r3))
192  return TCL_ERROR;
193  if (!get_3D_vector(interp,argv[4],r4))
194  return TCL_ERROR;
195 
196  r12 = r1 - r2;
197  r23 = r2 - r3;
198  r34 = r3 - r4;
199 
200  // Calculate the cross products and distances
201  Vector A = cross(r12,r23);
202  BigReal rA = A.length();
203  Vector B = cross(r23,r34);
204  BigReal rB = B.length();
205  Vector C = cross(r23,A);
206  BigReal rC = C.length();
207 
208  // Calculate the sin and cos
209  BigReal cos_phi = (A*B)/(rA*rB);
210  BigReal sin_phi = (C*B)/(rC*rB);
211 
212  Force f1,f2,f3;
213 
214  // Normalize B
215  rB = 1.0/rB;
216  B *= rB;
217 
218  // We first need to figure out whether the
219  // sin or cos form will be more stable. For this,
220  // just look at the value of phi
221  if (fabs(sin_phi) > 0.1)
222  {
223  // use the sin version to avoid 1/cos terms
224 
225  // Normalize A
226  rA = 1.0/rA;
227  A *= rA;
228  Vector dcosdA = rA*(cos_phi*A-B);
229  Vector dcosdB = rB*(cos_phi*B-A);
230 
231  K1 = -1/sin_phi;
232 
233  f1.x = K1*(r23.y*dcosdA.z - r23.z*dcosdA.y);
234  f1.y = K1*(r23.z*dcosdA.x - r23.x*dcosdA.z);
235  f1.z = K1*(r23.x*dcosdA.y - r23.y*dcosdA.x);
236 
237  f3.x = K1*(r23.z*dcosdB.y - r23.y*dcosdB.z);
238  f3.y = K1*(r23.x*dcosdB.z - r23.z*dcosdB.x);
239  f3.z = K1*(r23.y*dcosdB.x - r23.x*dcosdB.y);
240 
241  f2.x = K1*(r12.z*dcosdA.y - r12.y*dcosdA.z
242  + r34.y*dcosdB.z - r34.z*dcosdB.y);
243  f2.y = K1*(r12.x*dcosdA.z - r12.z*dcosdA.x
244  + r34.z*dcosdB.x - r34.x*dcosdB.z);
245  f2.z = K1*(r12.y*dcosdA.x - r12.x*dcosdA.y
246  + r34.x*dcosdB.y - r34.y*dcosdB.x);
247  }
248  else
249  {
250  // This angle is closer to 0 or 180 than it is to
251  // 90, so use the cos version to avoid 1/sin terms
252 
253  // Normalize C
254  rC = 1.0/rC;
255  C *= rC;
256  Vector dsindC = rC*(sin_phi*C-B);
257  Vector dsindB = rB*(sin_phi*B-C);
258 
259  K1 = 1/cos_phi;
260 
261  f1.x = K1*((r23.y*r23.y + r23.z*r23.z)*dsindC.x
262  - r23.x*r23.y*dsindC.y
263  - r23.x*r23.z*dsindC.z);
264  f1.y = K1*((r23.z*r23.z + r23.x*r23.x)*dsindC.y
265  - r23.y*r23.z*dsindC.z
266  - r23.y*r23.x*dsindC.x);
267  f1.z = K1*((r23.x*r23.x + r23.y*r23.y)*dsindC.z
268  - r23.z*r23.x*dsindC.x
269  - r23.z*r23.y*dsindC.y);
270 
271  f3 = cross(K1,dsindB,r23);
272 
273  f2.x = K1*(-(r23.y*r12.y + r23.z*r12.z)*dsindC.x
274  +(2.0*r23.x*r12.y - r12.x*r23.y)*dsindC.y
275  +(2.0*r23.x*r12.z - r12.x*r23.z)*dsindC.z
276  +dsindB.z*r34.y - dsindB.y*r34.z);
277  f2.y = K1*(-(r23.z*r12.z + r23.x*r12.x)*dsindC.y
278  +(2.0*r23.y*r12.z - r12.y*r23.z)*dsindC.z
279  +(2.0*r23.y*r12.x - r12.y*r23.x)*dsindC.x
280  +dsindB.x*r34.z - dsindB.z*r34.x);
281  f2.z = K1*(-(r23.x*r12.x + r23.y*r12.y)*dsindC.z
282  +(2.0*r23.z*r12.x - r12.z*r23.x)*dsindC.x
283  +(2.0*r23.z*r12.y - r12.z*r23.y)*dsindC.y
284  +dsindB.y*r34.x - dsindB.x*r34.y);
285  }
286 
287  Tcl_Obj* forcev[4];
288 
289  forcev[0] = obj_3D_vector(f1);
290  forcev[1] = obj_3D_vector(f2-f1);
291  forcev[2] = obj_3D_vector(f3-f2);
292  forcev[3] = obj_3D_vector(-f3);
293 
294  Tcl_SetObjResult(interp, Tcl_NewListObj(4, forcev));
295 
296  return TCL_OK;
297 }
static int get_3D_vector(Tcl_Interp *interp, Tcl_Obj *const list, Vector &result)
Definition: TclCommands.C:25
const BigReal A
Definition: Vector.h:64
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
BigReal z
Definition: Vector.h:66
BigReal length(void) const
Definition: Vector.h:169
BigReal x
Definition: Vector.h:66
static Tcl_Obj * obj_3D_vector(const Vector &v)
Definition: TclCommands.C:44
BigReal y
Definition: Vector.h:66
const BigReal B
double BigReal
Definition: common.h:114
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().

76 {
77  Vector r1, r2, r3, r12, r32;
78 
79  if (argc != 4)
80  return TCL_ERROR;
81  if (!get_3D_vector(interp,argv[1],r1))
82  return TCL_ERROR;
83  if (!get_3D_vector(interp,argv[2],r2))
84  return TCL_ERROR;
85  if (!get_3D_vector(interp,argv[3],r3))
86  return TCL_ERROR;
87  r12 = r1 - r2;
88  r32 = r3 - r2;
89  Tcl_SetObjResult(interp, Tcl_NewDoubleObj(
90  acos((r12*r32)/(r12.length()*r32.length()))*180/PI ));
91  return TCL_OK;
92 }
static int get_3D_vector(Tcl_Interp *interp, Tcl_Obj *const list, Vector &result)
Definition: TclCommands.C:25
Definition: Vector.h:64
BigReal length(void) const
Definition: Vector.h:169
#define PI
Definition: common.h:83
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().

59 {
60  Vector r1, r2;
61 
62  if (argc != 3)
63  return TCL_ERROR;
64  if (!get_3D_vector(interp,argv[1],r1))
65  return TCL_ERROR;
66  if (!get_3D_vector(interp,argv[2],r2))
67  return TCL_ERROR;
68  Tcl_SetObjResult(interp, Tcl_NewDoubleObj((r2-r1).length()));
69  return TCL_OK;
70 }
static int get_3D_vector(Tcl_Interp *interp, Tcl_Obj *const list, Vector &result)
Definition: TclCommands.C:25
Definition: Vector.h:64
int proc_getdihedral ( ClientData  ,
Tcl_Interp *  interp,
int  argc,
Tcl_Obj *const  argv[] 
)

Definition at line 97 of file TclCommands.C.

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

Referenced by tcl_vector_math_init().

98 {
99  BigReal rA, rB, rC;
100  Vector r1, r2, r3, r4, r12, r23, r34, A, B, C;
101 
102  if (argc != 5)
103  return TCL_ERROR;
104  if (!get_3D_vector(interp,argv[1],r1))
105  return TCL_ERROR;
106  if (!get_3D_vector(interp,argv[2],r2))
107  return TCL_ERROR;
108  if (!get_3D_vector(interp,argv[3],r3))
109  return TCL_ERROR;
110  if (!get_3D_vector(interp,argv[4],r4))
111  return TCL_ERROR;
112  r12 = r1 - r2;
113  r23 = r2 - r3;
114  r34 = r3 - r4;
115  A = cross(r12,r23);
116  B = cross(r23,r34);
117  C = cross(r23,A);
118  rA = A.length();
119  rB = B.length();
120  rC = C.length();
121  Tcl_SetObjResult(interp, Tcl_NewDoubleObj(
122  -atan2((C*B)/(rC*rB),(A*B)/(rA*rB))*180/PI ));
123  return TCL_OK;
124 }
static int get_3D_vector(Tcl_Interp *interp, Tcl_Obj *const list, Vector &result)
Definition: TclCommands.C:25
const BigReal A
Definition: Vector.h:64
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
BigReal length(void) const
Definition: Vector.h:169
#define PI
Definition: common.h:83
const BigReal B
double BigReal
Definition: common.h:114
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().

299  {
300 
301  // first import from TclVec.C stolen from VMD
302  Vec_Init(interp);
303 
304  Tcl_CreateObjCommand(interp, "getbond", proc_getbond,
305  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
306  Tcl_CreateObjCommand(interp, "getangle", proc_getangle,
307  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
308  Tcl_CreateObjCommand(interp, "getdihedral", proc_getdihedral,
309  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
310  Tcl_CreateObjCommand(interp, "anglegrad", proc_anglegrad,
311  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
312  Tcl_CreateObjCommand(interp, "dihedralgrad", proc_dihedralgrad,
313  (ClientData) NULL, (Tcl_CmdDeleteProc *) NULL);
314 
315  return TCL_OK;
316 }
int Vec_Init(Tcl_Interp *interp)
Definition: TclVec.C:643
int proc_dihedralgrad(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
Definition: TclCommands.C:180
int proc_getdihedral(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
Definition: TclCommands.C:97
int proc_anglegrad(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
Definition: TclCommands.C:130
int proc_getangle(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
Definition: TclCommands.C:75
int proc_getbond(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
Definition: TclCommands.C:58