NAMD
TclCommands.C
Go to the documentation of this file.
1 
7 #include <stdlib.h>
8 #ifndef _NO_MALLOC_H
9 #include <malloc.h>
10 #endif
11 #include <errno.h>
12 #include "TclCommands.h"
13 #include "Vector.h"
14 #include "NamdTypes.h"
15 
16 #ifdef NAMD_TCL
17 
18 #include <tcl.h>
19 
20 #include "Matrix4.C"
21 #include "TclVec.C"
22 
23 
24 // Get a 3-D vector from a TCL list
25 static int get_3D_vector(Tcl_Interp *interp, Tcl_Obj *const list, Vector &result)
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 }
41 
42 
43 // Append a 3-D vector to the result string
44 static Tcl_Obj* obj_3D_vector(const Vector &v)
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 }
54 
55 
56 // Function: getbond coor1 coor2
57 // Returns: the length of the bond formed by the two atoms (i.e., the distance between them)
58 int proc_getbond(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
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 }
71 
72 
73 // Function: getangle coor1 coor2 coor3
74 // Returns: the angle formed by the three atoms
75 int proc_getangle(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
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 }
93 
94 
95 // Function: getdihedral coor1 coor2 coor3 coor4
96 // Returns: the dihedral formed by the four atoms
97 int proc_getdihedral(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
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 }
125 
126 
127 // Function: anglegrad coor1 coor2 coor3
128 // Returns: a list of gradients for each atom
129 // The code was basically copied from ComputeAngles.C
130 int proc_anglegrad(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
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 }
175 
176 
177 // Function: dihedralgrad coor1 coor2 coor3 coor4
178 // Returns: a list of gradients for each atom
179 // The code was basically copied from ComputeDihedrals.C
180 int proc_dihedralgrad(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
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 }
298 
299 int tcl_vector_math_init(Tcl_Interp *interp) {
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 }
317 
318 
319 #endif
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
int Vec_Init(Tcl_Interp *interp)
Definition: TclVec.C:643
__device__ __forceinline__ float3 cross(const float3 v1, const float3 v2)
BigReal z
Definition: Vector.h:66
BigReal length(void) const
Definition: Vector.h:169
#define PI
Definition: common.h:81
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
BigReal x
Definition: Vector.h:66
static Tcl_Obj * obj_3D_vector(const Vector &v)
Definition: TclCommands.C:44
int tcl_vector_math_init(Tcl_Interp *interp)
Definition: TclCommands.C:299
BigReal y
Definition: Vector.h:66
const BigReal B
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
double BigReal
Definition: common.h:112
int proc_getbond(ClientData, Tcl_Interp *interp, int argc, Tcl_Obj *const argv[])
Definition: TclCommands.C:58