-
Notifications
You must be signed in to change notification settings - Fork 1
/
Copy pathenergy.c
86 lines (73 loc) · 2.24 KB
/
energy.c
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
/***********************************************************************
* energy.c
*
* (C) 2010 Nicolas Quesada, [email protected]
*
* This program 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 program 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 General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the Free Software
* Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
**********************************************************************/
/*! \file energy.c
\brief See documentation of energy.h - this file contains source code
only, with function prototypes in energy.h
*/
#include"energy.h"
#include"lj_params.h"
double
interaction_energy (double *x, double *y, double *z, int i, int n)
{
int j;
double U = 0;
double xij, yij, zij, rijm2, rijm6;
for (j = 0; j < i; j++)
{
xij = x[j] - x[i];
yij = y[j] - y[i];
zij = z[j] - z[i];
rijm2 = lj_sigma2 (i, j) / (xij * xij + yij * yij + zij * zij);
rijm6 = rijm2 * rijm2 * rijm2;
U += lj_epsilon (i, j) * rijm6 * (rijm6 - 1);
}
for (j = i + 1; j < n; j++)
{
xij = x[j] - x[i];
yij = y[j] - y[i];
zij = z[j] - z[i];
rijm2 = lj_sigma2 (i, j) / (xij * xij + yij * yij + zij * zij);
rijm6 = rijm2 * rijm2 * rijm2;
U += lj_epsilon (i, j) * rijm6 * (rijm6 - 1);
}
U = 4 * U;
return U;
}
double
cluster_energy (double *x, double *y, double *z, int n)
{
int i, j;
double U = 0;
double xij, yij, zij, rijm2, rijm6;
for (i = 0; i < n; i++)
{
for (j = 0; j < i; j++)
{
xij = x[j] - x[i];
yij = y[j] - y[i];
zij = z[j] - z[i];
rijm2 = lj_sigma2 (i, j) / (xij * xij + yij * yij + zij * zij);
rijm6 = rijm2 * rijm2 * rijm2;
U += lj_epsilon (i, j) * rijm6 * (rijm6 - 1);
}
}
U = 4 * U;
return U;
}