-
Notifications
You must be signed in to change notification settings - Fork 0
/
CalculateEnergy.h
233 lines (195 loc) · 8.59 KB
/
CalculateEnergy.h
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
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
/*******************************************************************************
GPU OPTIMIZED MONTE CARLO (GOMC) 2.31
Copyright (C) 2018 GOMC Group
A copy of the GNU General Public License can be found in the COPYRIGHT.txt
along with this program, also can be found at <http://www.gnu.org/licenses/>.
********************************************************************************/
#ifndef CALCULATEENERGY_H
#define CALCULATEENERGY_H
#include "BasicTypes.h"
#include "EnergyTypes.h"
#include "EwaldCached.h"
#include "Ewald.h"
#include "NoEwald.h"
#include "CellList.h"
#include <vector>
//
// CalculateEnergy.h
// Energy Calculation functions for Monte Carlo simulation
// Calculates using const references to a particular Simulation's members
// Brock Jackman Sep. 2013
//
// Updated to use radial-based intermolecular pressure
// Jason Mick Feb. 2014
//
// Updated with Mohammad B. to use neighbor list. Fixed various
// preexisting flaws in neighbor list codebase to make production ready.
// Jason Mick Feb. 2015
//
class StaticVals;
class System;
class Forcefield;
class Molecules;
class MoleculeLookup;
class MoleculeKind;
class Coordinates;
class COM;
class XYZArray;
class BoxDimensions;
namespace cbmc
{
class TrialMol;
}
class CalculateEnergy
{
public:
CalculateEnergy(StaticVals & stat, System & sys);
void Init(System & sys);
//! Calculates total energy/virial of all boxes in the system
SystemPotential SystemTotal() ;
//! Calculates total energy/virial of a single box in the system
SystemPotential BoxInter(SystemPotential potential,
XYZArray const& coords,
XYZArray const& com,
BoxDimensions const& boxAxes,
const uint box) ;
//! Calculate force and virial for the box
Virial ForceCalc(const uint box);
//! Calculates intermolecule energy of all boxes in the system
//! @param potential Copy of current energy structure to append result to
//! @param coords Particle coordinates to evaluate for
//! @param com Molecular centers of mass of system under evaluation
//! @param boxAxes Box Dimenions to evaluate in
//! @return System potential assuming no molecule changes
SystemPotential SystemInter(SystemPotential potential,
XYZArray const& coords,
XYZArray const& com,
BoxDimensions const& boxAxes) ;
//! Calculates intermolecular energy (LJ and coulomb) of a molecule
//! were it at molCoords.
//! @param potential Copy of current energy structure to append result to
//! @param molCoords Molecule coordinates
//! @param molIndex Index of molecule.
//! @param box Index of box molecule is in.
//! @param newCOM (optional) If COM has changed for new coordinate,
//! allows for that to be considered.
void MoleculeInter(Intermolecular &inter_LJ, Intermolecular &inter_coulomb,
XYZArray const& molCoords, const uint molIndex,
const uint box) const;
//! Calculates Nonbonded intra energy (LJ and coulomb )for
//! candidate positions
//! @param energy Return array, must be pre-allocated to size n
//! @param trialMol Partially built trial molecule.
//! @param trialPos Contains exactly n potential particle positions
//! @param partIndex Index of particle within the molecule
//! @param box Index of box molecule is in
//! @param trials Number of trials ot loop over in position array. (cbmc)
void ParticleNonbonded(double* inter, const cbmc::TrialMol& trialMol,
XYZArray const& trialPos,
const uint partIndex,
const uint box,
const uint trials) const;
//! Calculates Nonbonded inter energy (LJ and coulomb)for
//! candidate positions
//! @param energy Output Array, at least the size of trialpos
//! @param trialPos Array of candidate positions
//! @param partIndex Index of the particle within the molecule
//! @param molIndex Index of molecule
//! @param box Index of box molecule is in
//! @param trials Number of trials ot loop over in position array. (cbmc)
void ParticleInter(double* en, double *real,
XYZArray const& trialPos,
const uint partIndex,
const uint molIndex,
const uint box,
const uint trials) const;
void ParticleInterRange(double* en, double *real,
XYZArray const& trialPos,
const uint partIndex,
const uint molIndex,
const uint box,
const uint start,
const uint end) const;
//! Calculates the change in the TC from adding numChange atoms of a kind
//! @param box Index of box under consideration
//! @param kind Kind of particle being added or removed
//! @param add If removed: false (sign=-1); if added: true (sign=+1)
Intermolecular MoleculeTailChange(const uint box,
const uint kind,
const bool add) const;
//! Calculates intramolecular energy of a full molecule
void MoleculeIntra(const uint molIndex, const uint box, double *bondEn) const;
//! Calculates Nonbonded 1_3 intramolecule energy of a full molecule
//for Martini forcefield
double IntraEnergy_1_3(const double distSq, const uint atom1,
const uint atom2, const uint molIndex) const;
//! Calculates Nonbonded 1_4 intramolecule energy of a full molecule
//for Martini forcefield
double IntraEnergy_1_4(const double distSq, const uint atom1,
const uint atom2, const uint molIndex) const;
private:
//! Calculates full TC energy for one box in current system
void EnergyCorrection(SystemPotential& pot, BoxDimensions const& boxAxes,
const uint box) const;
//! Calculates full TC virial for one box in current system
void ForceCorrection(Virial& virial, BoxDimensions const& boxAxes,
const uint box) const;
//! Calculates bond vectors of a full molecule, stores them in vecs
void BondVectors(XYZArray & vecs,
MoleculeKind const& molKind,
const uint molIndex,
const uint box) const;
//! Calculates bond stretch intramolecular energy of a full molecule
void MolBond(double & energy,
MoleculeKind const& molKind,
XYZArray const& vecs,
const uint box) const;
//! Calculates angular bend intramolecular energy of a full molecule
void MolAngle(double & energy,
MoleculeKind const& molKind,
XYZArray const& vecs,
const uint box) const;
//! Calculates dihedral torsion intramolecular energy of a full molecule
void MolDihedral(double & energy,
MoleculeKind const& molKind,
XYZArray const& vecs,
const uint box) const;
//! Calculates Nonbonded 1_N intramolecule energy of a full molecule
void MolNonbond(double & energy,
MoleculeKind const& molKind,
const uint molIndex,
const uint box) const;
//! Calculates Nonbonded 1_4 intramolecule energy of a full molecule
void MolNonbond_1_4(double & energy,
MoleculeKind const& molKind,
const uint molIndex,
const uint box) const;
//! Calculates Nonbonded 1_3 intramolecule energy of a full molecule
//for Martini forcefield
void MolNonbond_1_3(double & energy,
MoleculeKind const& molKind,
const uint molIndex,
const uint box) const;
//! For particles in main coordinates array determines if they belong
//! to same molecule, using internal arrays.
bool SameMolecule(const uint p1, const uint p2) const
{
//We dont calculate the energy between two atom of same molecule or
uint pair1 = particleMol[p1];
uint pair2 = particleMol[p2];
return (pair1 == pair2);
}
const Forcefield& forcefield;
const Molecules& mols;
const Coordinates& currentCoords;
const MoleculeLookup& molLookup;
const BoxDimensions& currentAxes;
const COM& currentCOM;
const Ewald *calcEwald;
bool electrostatic, ewald;
std::vector<int> particleKind;
std::vector<int> particleMol;
std::vector<double> particleCharge;
const CellList& cellList;
};
#endif /*ENERGY_H*/