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 2 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

#	Line 1 | Line 1
1	<	#include "SRI.hpp"
2	<	#include "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	<	/**********************************************************************
18	<	*
19	<	* initialize vectors
20	<	*
21	<	***********************************************************************/
22	<
23	<	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) :
52	>	atom1_(atom1), atom2_(atom2), atom3_(atom3), atom4_(atom4), torsionType_(tt) { }
53
54	<	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;
54	>	void Torsion::calcForce(RealType& angle, bool doParticlePot) {
55
56	<	double r_cr1_sqr; /* the length of r_cr1 squared */
57	<	double r_cr2_sqr; /* the length of r_cr2 squared */
58	<
59	<	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];
56	>	Vector3d pos1 = atom1_->getPos();
57	>	Vector3d pos2 = atom2_->getPos();
58	>	Vector3d pos3 = atom3_->getPos();
59	>	Vector3d pos4 = atom4_->getPos();
60
61	<	c_p_a->getPos( aR );
62	<	c_p_b->getPos( bR );
63	<	c_p_c->getPos( cR );
48	<	c_p_d->getPos( dR );
61	>	Vector3d r21 = pos1 - pos2;
62	>	Vector3d r32 = pos2 - pos3;
63	>	Vector3d r43 = pos3 - pos4;
64
65	<	r_ab.x = bR[0] - aR[0];
66	<	r_ab.y = bR[1] - aR[1];
67	<	r_ab.z = bR[2] - aR[2];
68	<	r_ab.length  = sqrt((r_ab.x * r_ab.x + r_ab.y * r_ab.y + r_ab.z * r_ab.z));
65	>	//  Calculate the cross products and distances
66	>	Vector3d A = cross(r21, r32);
67	>	RealType rA = A.length();
68	>	Vector3d B = cross(r32, r43);
69	>	RealType rB = B.length();
70
71	<	r_cb.x = bR[0] - cR[0];
72	<	r_cb.y = bR[1] - cR[1];
73	<	r_cb.z = bR[2] - cR[2];
74	<	r_cb.length = sqrt((r_cb.x * r_cb.x + r_cb.y * r_cb.y + r_cb.z * r_cb.z));
75	<
76	<	r_cd.x = dR[0] - cR[0];
77	<	r_cd.y = dR[1] - cR[1];
78	<	r_cd.z = dR[2] - cR[2];
79	<	r_cd.length = sqrt((r_cd.x * r_cd.x + r_cd.y * r_cd.y + r_cd.z * r_cd.z));
80	<
81	<	r_cr1.x = r_ab.y * r_cb.z - r_cb.y * r_ab.z;
82	<	r_cr1.y = r_ab.z * r_cb.x - r_cb.z * r_ab.x;
83	<	r_cr1.z = r_ab.x * r_cb.y - r_cb.x * r_ab.y;
84	<	r_cr1_x2 = r_cr1.x * r_cr1.x;
85	<	r_cr1_y2 = r_cr1.y * r_cr1.y;
86	<	r_cr1_z2 = r_cr1.z * r_cr1.z;
87	<	r_cr1_sqr = r_cr1_x2 + r_cr1_y2 + r_cr1_z2;
88	<	r_cr1.length = sqrt(r_cr1_sqr);
89	<
90	<	r_cr2.x = r_cb.y * r_cd.z - r_cd.y * r_cb.z;
91	<	r_cr2.y = r_cb.z * r_cd.x - r_cd.z * r_cb.x;
92	<	r_cr2.z = r_cb.x * r_cd.y - r_cd.x * r_cb.y;
93	<	r_cr2_x2 = r_cr2.x * r_cr2.x;
94	<	r_cr2_y2 = r_cr2.y * r_cr2.y;
95	<	r_cr2_z2 = r_cr2.z * r_cr2.z;
96	<	r_cr2_sqr = r_cr2_x2 + r_cr2_y2 + r_cr2_z2;
97	<	r_cr2.length = sqrt(r_cr2_sqr);
98	<
99	<	r_cr1_r_cr2 = r_cr1.length * r_cr2.length;
100	<
101	<	/**********************************************************************
102	<	*
103	<	* dot product and angle calculations
104	<	*
105	<	***********************************************************************/
106	<
107	<	double cr1_dot_cr2; /* the dot product of the cr1 and cr2 vectors */
108	<	double cos_phi; /* the cosine of the torsion angle */
109	<
110	<	cr1_dot_cr2 = r_cr1.x * r_cr2.x + r_cr1.y * r_cr2.y + r_cr1.z * r_cr2.z;
111	<
112	<	cos_phi = cr1_dot_cr2 / r_cr1_r_cr2;
113	<
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);
71	>	/*
72	>	If either of the two cross product vectors is tiny, that means
73	>	the three atoms involved are colinear, and the torsion angle is
74	>	going to be undefined.  The easiest check for this problem is
75	>	to use the product of the two lengths.
76	>	*/
77	>	if (rA * rB < OpenMD::epsilon) return;
78	>
79	>	A.normalize();
80	>	B.normalize();
81	>
82	>	//  Calculate the sin and cos
83	>	RealType cos_phi = dot(A, B) ;
84	>	if (cos_phi > 1.0) cos_phi = 1.0;
85	>	if (cos_phi < -1.0) cos_phi = -1.0;
86	>
87	>	RealType dVdcosPhi;
88	>	torsionType_->calcForce(cos_phi, potential_, dVdcosPhi);
89	>	Vector3d f1 ;
90	>	Vector3d f2 ;
91	>	Vector3d f3 ;
92	>
93	>	Vector3d dcosdA = (cos_phi * A - B) /rA;
94	>	Vector3d dcosdB = (cos_phi * B - A) /rB;
95	>
96	>	f1 = dVdcosPhi * cross(r32, dcosdA);
97	>	f2 = dVdcosPhi * ( cross(r43, dcosdB) - cross(r21, dcosdA));
98	>	f3 = dVdcosPhi * cross(dcosdB, r32);
99	>
100	>	atom1_->addFrc(f1);
101	>	atom2_->addFrc(f2 - f1);
102	>	atom3_->addFrc(f3 - f2);
103	>	atom4_->addFrc(-f3);
104	>
105	>	if (doParticlePot) {
106	>	atom1_->addParticlePot(potential_);
107	>	atom2_->addParticlePot(potential_);
108	>	atom3_->addParticlePot(potential_);
109	>	atom4_->addParticlePot(potential_);
110	>	}
111	>
112	>	angle = acos(cos_phi) /M_PI * 180.0;
113	>	}
114		}

Comparing trunk/src/primitives/Torsion.cpp (property svn:keywords):
Revision 2 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
Revision 1879 by gezelter, Sun Jun 16 15:15:42 2013 UTC

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

Diff Legend

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