ViewVC Help

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

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

Comparing trunk/src/primitives/GhostBend.cpp (file contents):
Revision 2 by gezelter, Fri Sep 24 04:16:43 2004 UTC vs.
Revision 1953 by gezelter, Thu Dec 5 18:19:26 2013 UTC

#	Line 1 | Line 1
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	+	#include "primitives/GhostBend.hpp"
46	+	#include "primitives/DirectionalAtom.hpp"
47	+	namespace OpenMD {
48
49	<	#include <math.h>
50	<	#include <iostream>
51	<	#include <stdlib.h>
49	>	/**@todo still a lot left to improve*/
50	>	void GhostBend::calcForce(RealType& angle, bool doParticlePot) {
51	>	DirectionalAtom* ghostAtom = static_cast<DirectionalAtom*>(atoms_[1]);
52	>
53	>	Vector3d pos1 = atoms_[0]->getPos();
54	>	Vector3d pos2 = ghostAtom->getPos();
55
56	<	#include "simError.h"
57	<	#include "SRI.hpp"
58	<	#include "Atom.hpp"
56	>	Vector3d r21 = pos1 - pos2;
57	>	RealType d21 = r21.length();
58	>
59	>	RealType d21inv = 1.0 / d21;
60	>
61	>	// we need the transpose of A to get the lab fixed vector:
62	>	Vector3d r23 = ghostAtom->getA().transpose().getColumn(2);
63	>	RealType d23 = r23.length();
64	>
65	>	RealType d23inv = 1.0 / d23;
66	>
67	>	RealType cosTheta = dot(r21, r23) / (d21 * d23);
68
69	+	//check roundoff
70	+	if (cosTheta > 1.0) {
71	+	cosTheta = 1.0;
72	+	} else if (cosTheta < -1.0) {
73	+	cosTheta = -1.0;
74	+	}
75	+
76	+	RealType theta = acos(cosTheta);
77
78	<
12	<	GhostBend::GhostBend( Atom &a, Atom &b ){
13	<
14	<	c_p_a = &a;
15	<
16	<	if( !b.isDirectional() ){
78	>	RealType dVdTheta;
79
80	<	// if atom b is not directional, then bad things will happen
80	>	bendType_->calcForce(theta, potential_, dVdTheta);
81
82	<	sprintf( painCave.errMsg,
83	<	" Ghost Bend error: Atom # %d of type \"%s\" is not "
84	<	"directional.\n",
85	<	b.getIndex(),
86	<	b.getType() );
87	<	painCave.isFatal = 1;
88	<	simError();
89	<	}
82	>	RealType sinTheta = sqrt(1.0 - cosTheta * cosTheta);
83	>
84	>	if (fabs(sinTheta) < 1.0E-6) {
85	>	sinTheta = 1.0E-6;
86	>	}
87	>
88	>	RealType commonFactor1 = dVdTheta / sinTheta * d21inv;
89	>	RealType commonFactor2 = dVdTheta / sinTheta * d23inv;
90	>
91	>	Vector3d force1 = commonFactor1 * (r23 * d23inv - r21*d21inv*cosTheta);
92	>	Vector3d force3 = commonFactor2 * (r21 * d21inv - r23*d23inv*cosTheta);
93
94	<	atomB = ( DirectionalAtom* ) &b;
30	<
31	<	c_potential_E = 0.0;
94	>	// Total force in current bend is zero
95
96	<	}
96	>	atoms_[0]->addFrc(force1);
97	>	ghostAtom->addFrc(-force1);
98
99	+	ghostAtom->addTrq( cross(r23, force3) );
100	+	if(doParticlePot) {
101	+	atoms_[0]->addParticlePot(potential_);
102	+	ghostAtom->addParticlePot(potential_);
103	+	}
104
105	<	void GhostBend::calc_forces(){
106	<
107	<	double dx,dy,dz,gx,gy,gz,dx2,dy2,dz2,gx2,gy2,gz2;
108	<	double rij2, rkj2, riji2, rkji2, dot, denom, cosang, angl;
40	<
41	<	double sina2, sinai;
105	>	angle = theta /M_PI * 180.0;
106	>
107	>	}
108	>	} //end namespace OpenMD
109
43	–	double comf2, comf3, comf4;
44	–	double dcsidx, dcsidy, dcsidz, dcskdx, dcskdy, dcskdz;
45	–	// double dcsjdx, dcsjdy, dcsjdz;
46	–	double dadxi, dadyi, dadzi;
47	–	double dadxk, dadyk, dadzk;//, dadxj, dadyj, dadzj;
48	–	double daxi, dayi, dazi, daxk, dayk, dazk, daxj, dayj, dazj;
49	–	double u[3];
50	–
51	–	double aR[3], bR[3];
52	–	double aF[3], bF[3], bTrq[3];
53	–
54	–	c_p_a->getPos( aR );
55	–	atomB->getPos( bR );
56	–
57	–
58	–	dx = aR[0] - bR[0];
59	–	dy = aR[1] - bR[1];
60	–	dz = aR[2] - bR[2];
61	–
62	–	atomB->getU(u);
63	–
64	–	gx = u[0];
65	–	gy = u[1];
66	–	gz = u[2];
67	–
68	–	dx2 = dx * dx;
69	–	dy2 = dy * dy;
70	–	dz2 = dz * dz;
71	–
72	–	gx2 = gx * gx;
73	–	gy2 = gy * gy;
74	–	gz2 = gz * gz;
75	–
76	–	rij2 = dx2 + dy2 + dz2;
77	–	rkj2 = gx2 + gy2 + gz2;
78	–
79	–	riji2 = 1.0 / rij2;
80	–	rkji2 = 1.0 / rkj2;
81	–
82	–	dot = dx * gx + dy * gy + dz * gz;
83	–	denom = sqrt((riji2 * rkji2));
84	–	cosang = dot * denom;
85	–
86	–	if(cosang > 1.0)cosang = 1.0;
87	–	if(cosang < -1.0) cosang = -1.0;
88	–
89	–	angl = acos(cosang);
90	–	angl = angl * 180.0 / M_PI;
91	–
92	–	sina2 = 1.0 - cosang*cosang;
93	–	if(fabs(sina2) < 1.0E-12 ) sina2 = 1.0E-12;
94	–	sinai = 1.0 / sqrt(sina2);
95	–
96	–	comf2 = cosang * riji2;
97	–	comf3 = cosang * rkji2;
98	–	comf4 = bend_force(angl);
99	–
100	–	dcsidx = gx*denom - comf2*dx;
101	–	dcsidy = gy*denom - comf2*dy;
102	–	dcsidz = gz*denom - comf2*dz;
103	–
104	–	dcskdx = dx*denom - comf3*gx;
105	–	dcskdy = dy*denom - comf3*gy;
106	–	dcskdz = dz*denom - comf3*gz;
107	–
108	–	//   dcsjdx = -dcsidx - dcskdx;
109	–	//   dcsjdy = -dcsidy - dcskdy;
110	–	//   dcsjdz = -dcsidz - dcskdz;
111	–
112	–	dadxi = -sinai*dcsidx;
113	–	dadyi = -sinai*dcsidy;
114	–	dadzi = -sinai*dcsidz;
115	–
116	–	dadxk = -sinai*dcskdx;
117	–	dadyk = -sinai*dcskdy;
118	–	dadzk = -sinai*dcskdz;
119	–
120	–	//   dadxj = -dadxi - dadxk;
121	–	//   dadyj = -dadyi - dadyk;
122	–	//   dadzj = -dadzi - dadzk;
123	–
124	–	daxi = comf4*dadxi;
125	–	dayi = comf4*dadyi;
126	–	dazi = comf4*dadzi;
127	–
128	–	daxk = comf4*dadxk;
129	–	dayk = comf4*dadyk;
130	–	dazk = comf4*dadzk;
131	–
132	–	daxj = -daxi - daxk;
133	–	dayj = -dayi - dayk;
134	–	dazj = -dazi - dazk;
135	–
136	–	aF[0] = daxi;
137	–	aF[1] = dayi;
138	–	aF[2] = dazi;
139	–
140	–	bF[0] = daxj + daxk;
141	–	bF[1] = dayj + dayk;
142	–	bF[2] = dazj + dazk;
143	–
144	–	bTrq[0] = gy*dazk - gz*dayk;
145	–	bTrq[1] = gz*daxk - gx*dazk;
146	–	bTrq[2] = gx*dayk - gy*daxk;
147	–
148	–
149	–	c_p_a->addFrc( aF );
150	–	atomB->addFrc( bF );
151	–	atomB->addTrq( bTrq );
152	–
153	–	return;
154	–	}
155	–
156	–	void GhostBend::setConstants( double the_c1, double the_c2, double the_c3,
157	–	double the_Th0 ){
158	–	c1 = the_c1;
159	–	c2 = the_c2;
160	–	c3 = the_c3;
161	–	theta0 = the_Th0;
162	–	}
163	–
164	–
165	–	double GhostBend::bend_force( double theta ){
166	–
167	–	double dt, dt2;
168	–	double force;
169	–
170	–	dt = ( theta - theta0 ) * M_PI / 180.0;
171	–	dt2 = dt * dt;
172	–
173	–	c_potential_E = ( c1 * dt2 ) + ( c2 * dt ) + c3;
174	–	force = -( ( 2.0 * c1 * dt ) + c2 );
175	–	return force;
176	–	}

Comparing trunk/src/primitives/GhostBend.cpp (property svn:keywords):
Revision 2 by gezelter, Fri Sep 24 04:16:43 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