1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
|
/*
Project: Adun
Copyright (C) 2005 Michael Johnston & Jordi Villa-Freixa
Author: Michael Johnston
This application is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public
License as published by the Free Software Foundation; either
version 2 of the License, or (at your option) any later version.
This application is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
Library General Public License for more details.
You should have received a copy of the GNU General Public
License along with this library; if not, write to the Free
Software Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111 USA.
*/
/**
\note Using double quotes here with icc causes
many "function x referenced but not defined"
for external functions. I cannot figure this
out at all..
*/
#include <Base/AdForceFieldFunctions.h>
bool __HarmonicBondEnergyDebug__ = false;
bool __HarmonicBondForceDebug__ = false;
inline void AdHarmonicBondEnergy(double* bond, double **coordinates, double* bnd_pot)
{
register int i;
int atom_one, atom_two;
double eq_sep, bnd_cnst, forceMag;
Vector3D seperation_s;
//translate the information in bond
atom_one = (int)bond[0];
atom_two = (int)bond[1];
bnd_cnst = bond[2];
eq_sep = bond[3];
//calculate interatomic distance
for(i=3; --i>=0;)
*(seperation_s.vector + i) = coordinates[atom_two][i] - coordinates[atom_one][i];
//calculate the length of the seperation vector
Ad3DVectorLength(&seperation_s);
//calcualte acceleration magnitude
forceMag = -1*bnd_cnst*(seperation_s.length - eq_sep);
//calculate potential
*bnd_pot = *bnd_pot - forceMag*(seperation_s.length - eq_sep)*0.5;
#ifdef BASE_BONDED_DEBUG
if(__HarmonicBondEnergyDebug__)
{
fprintf(stderr, "%-6d%-6d%-12.5lf%-12.5lf%-12.5lf%-12.5lf\n",
atom_one,
atom_two,
bnd_cnst,
eq_sep,
seperation_s.length,
*bnd_pot);
}
#endif
}
inline void AdHarmonicBondForce(double* bond, double **coordinates, double **forces, double* bnd_pot)
{
register int i;
int atom_one, atom_two;
double forceMag, eq_sep, bnd_cnst, holder, rlength;
Vector3D seperation_s;
//translate the information in bond
atom_one = (int)bond[0];
atom_two = (int)bond[1];
bnd_cnst = bond[2];
eq_sep = bond[3];
//calculate interatomic distance
for(i=3; --i>=0;)
*(seperation_s.vector + i) = coordinates[atom_two][i] - coordinates[atom_one][i];
//calculate the length of the seperation vector
Ad3DVectorLength(&seperation_s);
//calcualte acceleration magnitude
forceMag = -1*bnd_cnst*(seperation_s.length - eq_sep);
//calculate potential
*bnd_pot = *bnd_pot - forceMag*(seperation_s.length - eq_sep)*0.5;
rlength = 1/seperation_s.length;
//calculate acceleration on atom one
for(i=0; i<3; i++)
{
holder = forceMag*seperation_s.vector[i]*rlength;
forces[atom_two][i] += holder;
forces[atom_one][i] -= holder;
}
#ifdef BASE_BONDED_DEBUG
if(__HarmonicBondForceDebug__)
{
fprintf(stderr, "%-6d%-6d%-12.5lf%-12.5lf%-12.5lf%-12.5lf%-12.5lf\n",
atom_one,
atom_two,
bnd_cnst,
eq_sep,
seperation_s.length,
*bnd_pot,
forceMag);
}
#endif
}
|
