-
Notifications
You must be signed in to change notification settings - Fork 0
/
Rho_HO.cpp
122 lines (107 loc) · 4.07 KB
/
Rho_HO.cpp
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
/*
* RhoHO.cpp
*
* Created on: May 17, 2012
* Author: dstuck
*/
#include "Rho_HO.h"
Rho_HO::Rho_HO(vector<double> w, int nFrozModes, int hFrozModes) : omega(w), lowFrozModes(nFrozModes), highFrozModes(hFrozModes) {
}
Rho_HO::Rho_HO(double w, int N) {
omega.clear();
for(int j=0; j<N; j++){
omega.push_back(w);
}
}
Rho_HO::~Rho_HO() {
// TODO Auto-generated destructor stub
}
double Rho_HO::GetRho(vector<Particle> slice1, vector<Particle> slice2, double eps) {
// This function is not called if using Levy Flights
double rho = 0;
vector<double> gamma;
if(omega.size() != slice1.size()) {
cout << "Error! Number of particles not equal to number of frequencies! Exiting with -99999." << endl;
return -99999;
}
for(int j=0; j<(int)slice1.size(); j++) {
for(int k=0; k<(int)slice1[j].pos.size(); k++) {
// rho += slice1[j].mass*omega[j]/2.0/eps * ((slice1[j].pos[k] - slice2[j].pos[k])*(slice1[j].pos[k] - slice2[j].pos[k])/eps/omega[j]/sinh(omega[j]*eps) + tanh(eps*omega[j]/2.0)/omega[j]/eps*(slice1[j].pos[k]*slice1[j].pos[k]+slice2[j].pos[k]*slice2[j].pos[k]));
rho += slice1[j].mass*omega[j]/2.0/eps * ((slice1[j].pos[k] - slice2[j].pos[k])*(slice1[j].pos[k] - slice2[j].pos[k])/sinh(omega[j]*eps) + tanh(eps*omega[j]/2.0)*(slice1[j].pos[k]*slice1[j].pos[k]+slice2[j].pos[k]*slice2[j].pos[k]));
}
}
return rho;
}
double Rho_HO::ModifyPotential(vector<Particle> part) {
double dV = 0.0;
int dim = part[0].pos.size();
for(int j=0; j<(int)part.size(); j++) {
for(int k=0; k<dim; k++) {
dV -= (part[j].pos[k])*(part[j].pos[k])/2.0*omega[j]*omega[j]*part[j].mass;
}
}
// cout << "\tV_HO =\t" << dV << endl;
// cout << dV << endl;
return dV;
}
double Rho_HO::Estimate(vector<Particle> slice1, vector<Particle> slice2, double eps, int P) {
/*
* Estimator taken from Whitfield and Martyna, J. Chem. Phys. 126, 074104 (2007) with missing factor of 1/2 added to term 2
*/
if(omega.size() != slice1.size()) {
cout << omega.size() << " frequencies for " << slice1.size() << " modes." << endl;
cout << "Error! Number of particles not equal to number of frequencies! Exiting with -99999." << endl;
return -99999;
}
double est = 0;
for(int j=0; j< lowFrozModes; j++) {
for(int k=0; k<(int)slice1[j].pos.size(); k++) {
est += omega[j]/2.0/tanh(eps*(double)P*omega[j]);
}
}
for(int j=(int)slice1.size()-highFrozModes; j<(int)slice1.size(); j++) {
for(int k=0; k<(int)slice1[j].pos.size(); k++) {
est += omega[j]/2.0/tanh(eps*(double)P*omega[j]);
}
}
for(int j=lowFrozModes; j<(int)slice1.size()-highFrozModes; j++) {
for(int k=0; k<(int)slice1[j].pos.size(); k++) {
est += omega[j]/2.0/tanh(eps*omega[j]);
est += -slice1[j].mass*omega[j]*omega[j]/2.0/tanh(eps*omega[j])/sinh(eps*omega[j])*(slice1[j].pos[k]-slice2[j].pos[k])*(slice1[j].pos[k]-slice2[j].pos[k]);
est += slice1[j].mass*omega[j]*omega[j]/2.0/cosh(eps*omega[j]/2.0)/cosh(eps*omega[j]/2.0)*(slice1[j].pos[k]*slice1[j].pos[k]);
}
}
est /= (double)P;
return est;
}
vector<double> Rho_HO::GetSpringLength(vector<Particle> part, double eps) {
cout << "Should not be calling this function!" << endl;
vector<double> r0;
for(int j=0; j<(int)part.size(); j++) {
r0.push_back(sqrt(tanh(eps*omega[j])/part[j].mass/omega[j]));
// r0.push_back(sqrt(sinh(eps*omega[j])/part[j].mass/omega[j]));
}
return r0;
}
vector<double> Rho_HO::GetLevyMean(Particle partI, Particle partF, double delta, double eps, int partNum) {
vector<double> mean;
double tDelI = tanh(eps*omega[partNum]);
double tDelF = tanh(delta*eps*omega[partNum]);
double sDelI = sinh(eps*omega[partNum]);
double sDelF = sinh(delta*eps*omega[partNum]);
for(int k = 0; k<(int)partI.pos.size(); k++) {
mean.push_back( (partI.pos[k]/sDelI + partF.pos[k]/sDelF) * tDelI*tDelF/(tDelI+tDelF) );
}
return mean;
}
double Rho_HO::GetLevySigma(double delta, double eps, int partNum) {
double sigma;
double tDelI = tanh(eps*omega[partNum]);
double tDelF = tanh(delta*eps*omega[partNum]);
sigma = sqrt(tDelI*tDelF/(tDelI+tDelF) / omega[partNum]);
return sigma;
}
string Rho_HO::GetType() {
string name = "HO";
return name;
}