ViewVC Help

View File | Revision Log | Show Annotations | View Changeset | Root Listing

root/OpenMD/trunk/src/primitives/Torsion.cpp

Comparing trunk/src/primitives/Torsion.cpp (file contents):
Revision 3 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
Revision 1953 by gezelter, Thu Dec 5 18:19:26 2013 UTC

#	Line 1 | Line 1
1	<	#include "primitives/SRI.hpp"
2	<	#include "primitives/Atom.hpp"
3	<	#include <math.h>
4	<	#include <iostream>
5	<	#include <stdlib.h>
1	>	/*
2	>	* Copyright (c) 2005 The University of Notre Dame. All Rights Reserved.
3	>	*
4	>	* The University of Notre Dame grants you ("Licensee") a
5	>	* non-exclusive, royalty free, license to use, modify and
6	>	* redistribute this software in source and binary code form, provided
7	>	* that the following conditions are met:
8	>	*
9	>	* 1. Redistributions of source code must retain the above copyright
10	>	*    notice, this list of conditions and the following disclaimer.
11	>	*
12	>	* 2. Redistributions in binary form must reproduce the above copyright
13	>	*    notice, this list of conditions and the following disclaimer in the
14	>	*    documentation and/or other materials provided with the
15	>	*    distribution.
16	>	*
17	>	* This software is provided "AS IS," without a warranty of any
18	>	* kind. All express or implied conditions, representations and
19	>	* warranties, including any implied warranty of merchantability,
20	>	* fitness for a particular purpose or non-infringement, are hereby
21	>	* excluded.  The University of Notre Dame and its licensors shall not
22	>	* be liable for any damages suffered by licensee as a result of
23	>	* using, modifying or distributing the software or its
24	>	* derivatives. In no event will the University of Notre Dame or its
25	>	* licensors be liable for any lost revenue, profit or data, or for
26	>	* direct, indirect, special, consequential, incidental or punitive
27	>	* damages, however caused and regardless of the theory of liability,
28	>	* arising out of the use of or inability to use software, even if the
29	>	* University of Notre Dame has been advised of the possibility of
30	>	* such damages.
31	>	*
32	>	* SUPPORT OPEN SCIENCE!  If you use OpenMD or its source code in your
33	>	* research, please cite the appropriate papers when you publish your
34	>	* work.  Good starting points are:
35	>	*
36	>	* [1]  Meineke, et al., J. Comp. Chem. 26, 252-271 (2005).
37	>	* [2]  Fennell & Gezelter, J. Chem. Phys. 124, 234104 (2006).
38	>	* [3]  Sun, Lin & Gezelter, J. Chem. Phys. 128, 234107 (2008).
39	>	* [4]  Kuang & Gezelter,  J. Chem. Phys. 133, 164101 (2010).
40	>	* [5]  Vardeman, Stocker & Gezelter, J. Chem. Theory Comput. 7, 834 (2011).
41	>	*/
42	>
43	>	#include "config.h"
44	>	#include <cmath>
45
46	<	void Torsion::set_atoms( Atom &a, Atom &b, Atom &c, Atom &d){
8	<	c_p_a = &a;
9	<	c_p_b = &b;
10	<	c_p_c = &c;
11	<	c_p_d = &d;
12	<	}
46	>	#include "primitives/Torsion.hpp"
47
48	+	namespace OpenMD {
49
50	<	void Torsion::calc_forces(){
51	<
52	<	/**********************************************************************
53	<	*
54	<	* initialize vectors
55	<	*
56	<	***********************************************************************/
57	<
58	<	vect r_ab; /* the vector whose origin is a and end is b */
24	<	vect r_cb; /* the vector whose origin is c and end is b */
25	<	vect r_cd; /* the vector whose origin is c and end is b */
26	<	vect r_cr1; /* the cross product of r_ab and r_cb */
27	<	vect r_cr2; /* the cross product of r_cb and r_cd */
50	>	Torsion::Torsion(Atom *atom1, Atom *atom2, Atom *atom3, Atom *atom4,
51	>	TorsionType *tt) : ShortRangeInteraction(),
52	>	torsionType_(tt) {
53	>	atoms_.resize(4);
54	>	atoms_[0] = atom1;
55	>	atoms_[1] = atom2;
56	>	atoms_[2] = atom3;
57	>	atoms_[3] = atom4;
58	>	}
59
60	<	double r_cr1_x2; /* the components of r_cr1 squared */
30	<	double r_cr1_y2;
31	<	double r_cr1_z2;
32	<
33	<	double r_cr2_x2; /* the components of r_cr2 squared */
34	<	double r_cr2_y2;
35	<	double r_cr2_z2;
60	>	void Torsion::calcForce(RealType& angle, bool doParticlePot) {
61
62	<	double r_cr1_sqr; /* the length of r_cr1 squared */
63	<	double r_cr2_sqr; /* the length of r_cr2 squared */
64	<
65	<	double r_cr1_r_cr2; /* the length of r_cr1 * length of r_cr2 */
41	<
42	<	double aR[3], bR[3], cR[3], dR[3];
43	<	double aF[3], bF[3], cF[3], dF[3];
62	>	Vector3d pos1 = atoms_[0]->getPos();
63	>	Vector3d pos2 = atoms_[1]->getPos();
64	>	Vector3d pos3 = atoms_[2]->getPos();
65	>	Vector3d pos4 = atoms_[3]->getPos();
66
67	<	c_p_a->getPos( aR );
68	<	c_p_b->getPos( bR );
69	<	c_p_c->getPos( cR );
48	<	c_p_d->getPos( dR );
67	>	Vector3d r21 = pos1 - pos2;
68	>	Vector3d r32 = pos2 - pos3;
69	>	Vector3d r43 = pos3 - pos4;
70
71	<	r_ab.x = bR[0] - aR[0];
72	<	r_ab.y = bR[1] - aR[1];
73	<	r_ab.z = bR[2] - aR[2];
74	<	r_ab.length  = sqrt((r_ab.x * r_ab.x + r_ab.y * r_ab.y + r_ab.z * r_ab.z));
71	>	//  Calculate the cross products and distances
72	>	Vector3d A = cross(r21, r32);
73	>	RealType rA = A.length();
74	>	Vector3d B = cross(r32, r43);
75	>	RealType rB = B.length();
76
77	<	r_cb.x = bR[0] - cR[0];
78	<	r_cb.y = bR[1] - cR[1];
79	<	r_cb.z = bR[2] - cR[2];
80	<	r_cb.length = sqrt((r_cb.x * r_cb.x + r_cb.y * r_cb.y + r_cb.z * r_cb.z));
81	<
82	<	r_cd.x = dR[0] - cR[0];
83	<	r_cd.y = dR[1] - cR[1];
84	<	r_cd.z = dR[2] - cR[2];
85	<	r_cd.length = sqrt((r_cd.x * r_cd.x + r_cd.y * r_cd.y + r_cd.z * r_cd.z));
86	<
87	<	r_cr1.x = r_ab.y * r_cb.z - r_cb.y * r_ab.z;
88	<	r_cr1.y = r_ab.z * r_cb.x - r_cb.z * r_ab.x;
89	<	r_cr1.z = r_ab.x * r_cb.y - r_cb.x * r_ab.y;
90	<	r_cr1_x2 = r_cr1.x * r_cr1.x;
91	<	r_cr1_y2 = r_cr1.y * r_cr1.y;
92	<	r_cr1_z2 = r_cr1.z * r_cr1.z;
93	<	r_cr1_sqr = r_cr1_x2 + r_cr1_y2 + r_cr1_z2;
94	<	r_cr1.length = sqrt(r_cr1_sqr);
95	<
96	<	r_cr2.x = r_cb.y * r_cd.z - r_cd.y * r_cb.z;
97	<	r_cr2.y = r_cb.z * r_cd.x - r_cd.z * r_cb.x;
98	<	r_cr2.z = r_cb.x * r_cd.y - r_cd.x * r_cb.y;
99	<	r_cr2_x2 = r_cr2.x * r_cr2.x;
100	<	r_cr2_y2 = r_cr2.y * r_cr2.y;
101	<	r_cr2_z2 = r_cr2.z * r_cr2.z;
102	<	r_cr2_sqr = r_cr2_x2 + r_cr2_y2 + r_cr2_z2;
103	<	r_cr2.length = sqrt(r_cr2_sqr);
104	<
105	<	r_cr1_r_cr2 = r_cr1.length * r_cr2.length;
106	<
107	<	/**********************************************************************
108	<	*
109	<	* dot product and angle calculations
110	<	*
111	<	***********************************************************************/
112	<
113	<	double cr1_dot_cr2; /* the dot product of the cr1 and cr2 vectors */
114	<	double cos_phi; /* the cosine of the torsion angle */
115	<
116	<	cr1_dot_cr2 = r_cr1.x * r_cr2.x + r_cr1.y * r_cr2.y + r_cr1.z * r_cr2.z;
117	<
118	<	cos_phi = cr1_dot_cr2 / r_cr1_r_cr2;
119	<
98	<	/* adjust for the granularity of the numbers for angles near 0 or pi */
99	<
100	<	if(cos_phi > 1.0) cos_phi = 1.0;
101	<	if(cos_phi < -1.0) cos_phi = -1.0;
102	<
103	<
104	<	/********************************************************************
105	<	*
106	<	* This next section calculates derivatives needed for the force
107	<	* calculation
108	<	*
109	<	********************************************************************/
110	<
111	<
112	<	/* the derivatives of cos phi with respect to the x, y,
113	<	and z components of vectors cr1 and cr2. */
114	<	double d_cos_dx_cr1;
115	<	double d_cos_dy_cr1;
116	<	double d_cos_dz_cr1;
117	<	double d_cos_dx_cr2;
118	<	double d_cos_dy_cr2;
119	<	double d_cos_dz_cr2;
120	<
121	<	d_cos_dx_cr1 = r_cr2.x / r_cr1_r_cr2 - (cos_phi * r_cr1.x) / r_cr1_sqr;
122	<	d_cos_dy_cr1 = r_cr2.y / r_cr1_r_cr2 - (cos_phi * r_cr1.y) / r_cr1_sqr;
123	<	d_cos_dz_cr1 = r_cr2.z / r_cr1_r_cr2 - (cos_phi * r_cr1.z) / r_cr1_sqr;
124	<
125	<	d_cos_dx_cr2 = r_cr1.x / r_cr1_r_cr2 - (cos_phi * r_cr2.x) / r_cr2_sqr;
126	<	d_cos_dy_cr2 = r_cr1.y / r_cr1_r_cr2 - (cos_phi * r_cr2.y) / r_cr2_sqr;
127	<	d_cos_dz_cr2 = r_cr1.z / r_cr1_r_cr2 - (cos_phi * r_cr2.z) / r_cr2_sqr;
128	<
129	<	/***********************************************************************
130	<	*
131	<	* Calculate the actual forces and place them in the atoms.
132	<	*
133	<	***********************************************************************/
134	<
135	<	double force; /*the force scaling factor */
136	<
137	<	force = torsion_force(cos_phi);
138	<
139	<	aF[0] = force * (d_cos_dy_cr1 * r_cb.z - d_cos_dz_cr1 * r_cb.y);
140	<	aF[1] = force * (d_cos_dz_cr1 * r_cb.x - d_cos_dx_cr1 * r_cb.z);
141	<	aF[2] = force * (d_cos_dx_cr1 * r_cb.y - d_cos_dy_cr1 * r_cb.x);
142	<
143	<	bF[0] = force * (  d_cos_dy_cr1 * (r_ab.z - r_cb.z)
144	<	- d_cos_dy_cr2 *  r_cd.z
145	<	+ d_cos_dz_cr1 * (r_cb.y - r_ab.y)
146	<	+ d_cos_dz_cr2 *  r_cd.y);
147	<	bF[1] = force * (  d_cos_dx_cr1 * (r_cb.z - r_ab.z)
148	<	+ d_cos_dx_cr2 *  r_cd.z
149	<	+ d_cos_dz_cr1 * (r_ab.x - r_cb.x)
150	<	- d_cos_dz_cr2 *  r_cd.x);
151	<	bF[2] = force * (  d_cos_dx_cr1 * (r_ab.y - r_cb.y)
152	<	- d_cos_dx_cr2 *  r_cd.y
153	<	+ d_cos_dy_cr1 * (r_cb.x - r_ab.x)
154	<	+ d_cos_dy_cr2 *  r_cd.x);
155	<
156	<	cF[0] = force * (- d_cos_dy_cr1 *  r_ab.z
157	<	- d_cos_dy_cr2 * (r_cb.z - r_cd.z)
158	<	+ d_cos_dz_cr1 *  r_ab.y
159	<	- d_cos_dz_cr2 * (r_cd.y - r_cb.y));
160	<	cF[1] = force * (  d_cos_dx_cr1 *  r_ab.z
161	<	- d_cos_dx_cr2 * (r_cd.z - r_cb.z)
162	<	- d_cos_dz_cr1 *  r_ab.x
163	<	- d_cos_dz_cr2 * (r_cb.x - r_cd.x));
164	<	cF[2] = force * (- d_cos_dx_cr1 *  r_ab.y
165	<	- d_cos_dx_cr2 * (r_cb.y - r_cd.y)
166	<	+ d_cos_dy_cr1 *  r_ab.x
167	<	- d_cos_dy_cr2 * (r_cd.x - r_cb.x));
168	<
169	<	dF[0] = force * (d_cos_dy_cr2 * r_cb.z - d_cos_dz_cr2 * r_cb.y);
170	<	dF[1] = force * (d_cos_dz_cr2 * r_cb.x - d_cos_dx_cr2 * r_cb.z);
171	<	dF[2] = force * (d_cos_dx_cr2 * r_cb.y - d_cos_dy_cr2 * r_cb.x);
172	<
173	<
174	<	c_p_a->addFrc(aF);
175	<	c_p_b->addFrc(bF);
176	<	c_p_c->addFrc(cF);
177	<	c_p_d->addFrc(dF);
77	>	/*
78	>	If either of the two cross product vectors is tiny, that means
79	>	the three atoms involved are colinear, and the torsion angle is
80	>	going to be undefined.  The easiest check for this problem is
81	>	to use the product of the two lengths.
82	>	*/
83	>	if (rA * rB < OpenMD::epsilon) return;
84	>
85	>	A.normalize();
86	>	B.normalize();
87	>
88	>	//  Calculate the sin and cos
89	>	RealType cos_phi = dot(A, B) ;
90	>	if (cos_phi > 1.0) cos_phi = 1.0;
91	>	if (cos_phi < -1.0) cos_phi = -1.0;
92	>
93	>	RealType dVdcosPhi;
94	>	torsionType_->calcForce(cos_phi, potential_, dVdcosPhi);
95	>	Vector3d f1 ;
96	>	Vector3d f2 ;
97	>	Vector3d f3 ;
98	>
99	>	Vector3d dcosdA = (cos_phi * A - B) /rA;
100	>	Vector3d dcosdB = (cos_phi * B - A) /rB;
101	>
102	>	f1 = dVdcosPhi * cross(r32, dcosdA);
103	>	f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA));
104	>	f3 = dVdcosPhi * cross(dcosdB, r32);
105	>
106	>	atoms_[0]->addFrc(f1);
107	>	atoms_[1]->addFrc(f2 - f1);
108	>	atoms_[2]->addFrc(f3 - f2);
109	>	atoms_[3]->addFrc(-f3);
110	>
111	>	if (doParticlePot) {
112	>	atoms_[0]->addParticlePot(potential_);
113	>	atoms_[1]->addParticlePot(potential_);
114	>	atoms_[2]->addParticlePot(potential_);
115	>	atoms_[3]->addParticlePot(potential_);
116	>	}
117	>
118	>	angle = acos(cos_phi) /M_PI * 180.0;
119	>	}
120		}

Comparing trunk/src/primitives/Torsion.cpp (property svn:keywords):
Revision 3 by tim, Fri Sep 24 16:27:58 2004 UTC vs.
Revision 1953 by gezelter, Thu Dec 5 18:19:26 2013 UTC

#	Line 0 | Line 1
1	+	Author Id Revision Date

Diff Legend

–	Removed lines
+	Added lines
<	Changed lines
>	Changed lines