#include <ComputeAngles.h>
|
|
Definition at line 19 of file ComputeAngles.h. 00019 { size = 3 };
|
|
|
Definition at line 47 of file ComputeAngles.h. 00047 { angleEnergyIndex, TENSOR(virialIndex), reductionDataSize };
|
|
|
Definition at line 48 of file ComputeAngles.h. 00048 { reductionChecksumLabel = REDUCTION_ANGLE_CHECKSUM };
|
|
|
Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved. Definition at line 12 of file ComputeAngles.inl. 00012 { ; }
|
|
||||||||||||||||
|
Definition at line 14 of file ComputeAngles.inl. References atomID, TupleSignature::offset, TupleSignature::tupleParamType, and value. 00014 {
00015 atomID[0] = atom0;
00016 atomID[1] = atom0 + sig->offset[0];
00017 atomID[2] = atom0 + sig->offset[1];
00018 value = &v[sig->tupleParamType];
00019 }
|
|
||||||||||||
|
Definition at line 21 of file ComputeAngles.inl. References Angle, angle::angle_type, angle::atom1, angle::atom2, angle::atom3, atomID, and value. 00022 {
00023 atomID[0] = a->atom1;
00024 atomID[1] = a->atom2;
00025 atomID[2] = a->atom3;
00026 value = &v[a->angle_type];
00027 }
|
|
||||||||||||||||
|
Definition at line 29 of file ComputeAngles.inl. References atomID, and AtomID. 00030 {
00031 if (atom0 > atom2) { // Swap end atoms so lowest is first!
00032 AtomID tmp = atom2; atom2 = atom0; atom0 = tmp;
00033 }
00034 atomID[0] = atom0;
00035 atomID[1] = atom1;
00036 atomID[2] = atom2;
00037 }
|
|
|
Definition at line 55 of file ComputeAngles.h. 00055 { };
|
|
||||||||||||
|
Definition at line 72 of file ComputeAngles.C. References BigReal, DebugM, Lattice::delta(), TuplePatchElem::f, Force, AngleValue::k, AngleValue::k_ub, Patch::lattice, Vector::length(), localIndex, TuplePatchElem::p, p, CompAtom::partition, CompAtom::position, Position, pp_clamp(), pp_reduction(), pressureProfileMin, pressureProfileSlabs, pressureProfileThickness, AngleValue::r_ub, Vector::rlength(), AngleValue::theta0, value, Vector::x, TuplePatchElem::x, Vector::y, and Vector::z. 00073 {
00074 DebugM(3, "::computeForce() localIndex = " << localIndex[0] << " "
00075 << localIndex[1] << " " << localIndex[2] << std::endl);
00076
00077 const Position & pos1 = p[0]->x[localIndex[0]].position;
00078 const Lattice & lattice = p[0]->p->lattice;
00079 const Position & pos2 = p[1]->x[localIndex[1]].position;
00080 Vector r12 = lattice.delta(pos1,pos2);
00081 register BigReal d12inv = r12.rlength();
00082 const Position & pos3 = p[2]->x[localIndex[2]].position;
00083 Vector r32 = lattice.delta(pos3,pos2);
00084 register BigReal d32inv = r32.rlength();
00085
00086 BigReal cos_theta = (r12*r32)*(d12inv*d32inv);
00087 // Make sure that the cosine value is acceptable. With roundoff, you
00088 // can get values like 1.0+2e-16, which makes acos puke. So instead,
00089 // just set these kinds of values to exactly 1.0
00090 if (cos_theta > 1.0) cos_theta = 1.0;
00091 else if (cos_theta < -1.0) cos_theta = -1.0;
00092
00093 BigReal k = value->k * scale;
00094 BigReal theta0 = value->theta0;
00095
00096 // Get theta
00097 BigReal theta = acos(cos_theta);
00098
00099 // Compare it to the rest angle
00100 BigReal diff = theta - theta0;
00101
00102 // Add the energy from this angle to the total energy
00103 BigReal energy = k *diff*diff;
00104
00105 // Normalize vector r12 and r32
00106 //BigReal d12inv = 1. / d12;
00107 //BigReal d32inv = 1. / d32;
00108
00109 // Calculate constant factor 2k(theta-theta0)/sin(theta)
00110 BigReal sin_theta = sqrt(1.0 - cos_theta*cos_theta);
00111 if ( sin_theta < 1.e-6 ) {
00112 // catch case where bonds are parallel
00113 // this gets the force approximately right for theta0 of 0 or pi
00114 // and at least avoids small division for other values
00115 if ( diff < 0. ) diff = 2.0 * k;
00116 else diff = -2.0 * k;
00117 } else { diff *= (-2.0* k) / sin_theta; }
00118 BigReal c1 = diff * d12inv;
00119 BigReal c2 = diff * d32inv;
00120
00121 // Calculate the actual forces
00122 Force force1 = c1*(r12*(d12inv*cos_theta) - r32*d32inv);
00123 Force force2 = force1;
00124 Force force3 = c2*(r32*(d32inv*cos_theta) - r12*d12inv);
00125 force2 += force3; force2 *= -1;
00126
00127 // Check to see if we need to do the Urey-Bradley term
00128 if (value->k_ub)
00129 {
00130 // Non-zero k_ub value, so calculate the harmonic
00131 // potential between the 1-3 atoms
00132 BigReal k_ub = value->k_ub;
00133 BigReal r_ub = value->r_ub;
00134 Vector r13 = r12 - r32;
00135 BigReal d13 = r13.length();
00136 diff = d13- r_ub;
00137
00138 energy += k_ub *diff*diff;
00139
00140 diff *= -2.0*k_ub / d13;
00141 r13 *= diff;
00142
00143 force1 += r13;
00144 force3 -= r13;
00145 }
00146
00147 p[0]->f[localIndex[0]].x += force1.x;
00148 p[0]->f[localIndex[0]].y += force1.y;
00149 p[0]->f[localIndex[0]].z += force1.z;
00150
00151 p[1]->f[localIndex[1]].x += force2.x;
00152 p[1]->f[localIndex[1]].y += force2.y;
00153 p[1]->f[localIndex[1]].z += force2.z;
00154
00155 p[2]->f[localIndex[2]].x += force3.x;
00156 p[2]->f[localIndex[2]].y += force3.y;
00157 p[2]->f[localIndex[2]].z += force3.z;
00158
00159 DebugM(3, "::computeForce() -- ending with delta energy " << energy << std::endl);
00160
00161 reduction[angleEnergyIndex] += energy;
00162 reduction[virialIndex_XX] += ( force1.x * r12.x + force3.x * r32.x );
00163 reduction[virialIndex_XY] += ( force1.x * r12.y + force3.x * r32.y );
00164 reduction[virialIndex_XZ] += ( force1.x * r12.z + force3.x * r32.z );
00165 reduction[virialIndex_YX] += ( force1.y * r12.x + force3.y * r32.x );
00166 reduction[virialIndex_YY] += ( force1.y * r12.y + force3.y * r32.y );
00167 reduction[virialIndex_YZ] += ( force1.y * r12.z + force3.y * r32.z );
00168 reduction[virialIndex_ZX] += ( force1.z * r12.x + force3.z * r32.x );
00169 reduction[virialIndex_ZY] += ( force1.z * r12.y + force3.z * r32.y );
00170 reduction[virialIndex_ZZ] += ( force1.z * r12.z + force3.z * r32.z );
00171
00172 if (pressureProfileData) {
00173 BigReal z1 = p[0]->x[localIndex[0]].position.z;
00174 BigReal z2 = p[1]->x[localIndex[1]].position.z;
00175 BigReal z3 = p[2]->x[localIndex[2]].position.z;
00176 int n1 = (int)floor((z1-pressureProfileMin)/pressureProfileThickness);
00177 int n2 = (int)floor((z2-pressureProfileMin)/pressureProfileThickness);
00178 int n3 = (int)floor((z3-pressureProfileMin)/pressureProfileThickness);
00179 pp_clamp(n1, pressureProfileSlabs);
00180 pp_clamp(n2, pressureProfileSlabs);
00181 pp_clamp(n3, pressureProfileSlabs);
00182 int p1 = p[0]->x[localIndex[0]].partition;
00183 int p2 = p[1]->x[localIndex[1]].partition;
00184 int p3 = p[2]->x[localIndex[2]].partition;
00185 int pn = pressureProfileAtomTypes;
00186 pp_reduction(pressureProfileSlabs, n1, n2,
00187 p1, p2, pn,
00188 force1.x * r12.x, force1.y * r12.y, force1.z * r12.z,
00189 pressureProfileData);
00190 pp_reduction(pressureProfileSlabs, n3, n2,
00191 p3, p2, pn,
00192 force3.x * r32.x, force3.y * r32.y, force3.z * r32.z,
00193 pressureProfileData);
00194 }
00195 }
|
|
||||||||||||||||||||
|
Definition at line 57 of file ComputeAngles.C. References Angle, int32, and NAMD_die(). 00058 {
00059 #ifdef MEM_OPT_VERSION
00060 NAMD_die("Should not be called in AngleElem::getMoleculePointers in memory optimized version!");
00061 #else
00062 *count = mol->numAngles;
00063 *byatom = mol->anglesByAtom;
00064 *structarray = mol->angles;
00065 #endif
00066 }
|
|
||||||||||||
|
Definition at line 68 of file ComputeAngles.C. References Parameters::angle_array. 00068 {
00069 *v = p->angle_array;
00070 }
|
|
||||||||||||||||
|
Definition at line 29 of file ComputeAngles.h. References AtomSignature::angleCnt, and AtomSignature::angleSigs.
|
|
|
Definition at line 43 of file ComputeAngles.h.
|
|
||||||||||||||||
|
|
|
|
Definition at line 45 of file ComputeAngles.inl. References atomID. 00046 {
00047 return (atomID[0] < a.atomID[0] ||
00048 (atomID[0] == a.atomID[0] &&
00049 (atomID[1] < a.atomID[1] ||
00050 (atomID[1] == a.atomID[1] &&
00051 atomID[2] < a.atomID[2]) )));
00052 }
|
|
|
Definition at line 39 of file ComputeAngles.inl. References atomID. 00040 {
00041 return (a.atomID[0] == atomID[0] && a.atomID[1] == atomID[1] &&
00042 a.atomID[2] == atomID[2]);
00043 }
|
|
||||||||||||
|
Definition at line 198 of file ComputeAngles.C. References ADD_TENSOR, SubmitReduction::item(), REDUCTION_ANGLE_ENERGY, REDUCTION_VIRIAL_NORMAL, and virialIndex. 00199 {
00200 reduction->item(REDUCTION_ANGLE_ENERGY) += data[angleEnergyIndex];
00201 ADD_TENSOR(reduction,REDUCTION_VIRIAL_NORMAL,data,virialIndex);
00202 }
|
|
|
Definition at line 20 of file ComputeAngles.h. Referenced by AngleElem(), operator<(), and operator==(). |
|
|
Definition at line 21 of file ComputeAngles.h. Referenced by computeForce(). |
|
|
Definition at line 22 of file ComputeAngles.h. Referenced by computeForce(). |
|
|
Definition at line 52 of file ComputeAngles.C. |
|
|
Definition at line 54 of file ComputeAngles.C. Referenced by computeForce(). |
|
|
Copyright (c) 1995, 1996, 1997, 1998, 1999, 2000 by The Board of Trustees of the University of Illinois. All rights reserved. Definition at line 51 of file ComputeAngles.C. Referenced by computeForce(). |
|
|
Definition at line 53 of file ComputeAngles.C. Referenced by computeForce(). |
|
|
Definition at line 23 of file ComputeAngles.h. |
|
|
Definition at line 41 of file ComputeAngles.h. Referenced by AngleElem(), and computeForce(). |
1.3.9.1