Skip to content

Commit

Permalink
Added internals using modified library.
Browse files Browse the repository at this point in the history
  • Loading branch information
dstuck committed Mar 11, 2015
1 parent e567d4f commit 0880537
Show file tree
Hide file tree
Showing 73 changed files with 15,953 additions and 294 deletions.
921 changes: 830 additions & 91 deletions CoordUtil.cpp

Large diffs are not rendered by default.

27 changes: 24 additions & 3 deletions CoordUtil.h
Original file line number Diff line number Diff line change
Expand Up @@ -15,25 +15,34 @@
#include <vector>
#include <map>
#include "armadillo"
#include "optking/molecule.h"
using namespace std;

class CoordUtil {
public:
CoordUtil();
CoordUtil(int, int, vector< vector< vector<double> > >, vector<double>, vector<double>, vector< vector<double> >, vector<string>, vector<int>, vector< vector<int> >, string, string,bool,bool);
CoordUtil(int, int, bool, vector< vector< vector<double> > >, vector<double>, vector<double>, vector< vector<double> >, vector<string>, vector<int>, vector< vector<int> >, string, string,bool,bool);
virtual ~CoordUtil();

void initInternals();
vector< vector<double> > normalModeToCart(vector<Particle>);
vector< vector<double> > normalModeToCart(vector<double>); //Inefficient, for debug
vector<double> cartToNormalMode(vector< vector<double> >);
void MatchModes(CoordUtil*);
void MakeScalingVec(CoordUtil*);
void initMass();
int atomStrToInt(string);
CoordUtil* Clone();

bool useGuess;
bool internals;
bool readOmega;
bool readGeom;
bool scaleGeom;
int numModes;
int numPart;
int dim;
std::string tinkerName;
std::string outFileName;
std::string prmName;
map<string,double> atomToMass;
vector<string> atomType;
Expand All @@ -43,8 +52,20 @@ class CoordUtil {
vector<double> omega;
vector< vector<int> > connectivity;
vector< vector<double> > initCart;
vector< vector< vector<double> > > normModes;
vector< vector<double> > guessCarts;
vector<Particle> guessPart;
vector< vector< vector<double> > > normModes; //TODO: remove this in favor of armaModes
arma::vec scalingVec;
arma::vec atomMass;
arma::mat fragExtraB;
arma::mat internalModes;
arma::mat armaInvModes;
arma::mat sqrtInvMass; //TODO: Don't really need this as class data
//arma::mat equibBInv;
arma::mat equibBMat;
arma::mat delocIntCos;
//arma::mat wMat; // Weight matrix to normalize equibBInv, not sure if this is right...
opt::MOLECULE * workingMol;
};

#endif /* COORDUTIL_H_ */
132 changes: 132 additions & 0 deletions Est_AutoCorr.cpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,132 @@
/*
* Est_AutoCorr.cpp
*
* Created on: March 3, 2015
* Author: dstuck
*/

#include "Est_AutoCorr.h"
#include "armadillo"

Est_AutoCorr::Est_AutoCorr() {
}

Est_AutoCorr::Est_AutoCorr(int maxH) : maxHist(maxH) {
//hist = new deque<double>();
//totalStats = Stats();
Stats tempStats = Stats();
for(int i=0; i<maxH; i++) {
corrStats.push_back(tempStats);
}
nSamples=0;
}

Est_AutoCorr::~Est_AutoCorr() {
//delete hist;
}

void Est_AutoCorr::Reset() {
for(int i=0;i<corrStats.size();i++) {
corrStats[i].Reset();
}
//delete hist;
//hist = new deque<double>();
hist.clear();
nSamples = 0;

}

void Est_AutoCorr::AddVal(double newVal) {
//Add val to history
hist.push_front(newVal);
if(hist.size()>maxHist) {
hist.pop_back();
}
nSamples++;

totalStats.AddVal(newVal);
double newMean = totalStats.GetMean();
newVal -= newMean;
int count=0;
for(deque<double>::iterator iter=hist.begin(); iter!=hist.end(); iter++) {
corrStats[count].AddVal((*iter-newMean) * newVal);
count++;
}
}

vector<double> Est_AutoCorr::GetCorr() {
if(nSamples<maxHist) {
cout << "DES Warning: number of samples less than max history in Est_AutoCorr" << endl;
cout << nSamples << " < " << maxHist << endl;
}
vector<double> corrFunc (maxHist, 0.0);
//double var = totalStats.GetVariance();
double var = corrStats[0].GetMean()/double(nSamples);
//var /= min(maxHist,nSamples);
// cout << "totalStats.GetVariance() =" << totalStats.GetVariance() << endl;
// cout << "corrStats[0].GetVariance() = " << corrStats[0].GetVariance() << endl;
// cout << "var = " << var << endl;
// cout << "corrStats[0] mean = " << corrStats[0].GetMean() << endl;;
for(int i=0; i<min(maxHist,nSamples); i++) {
corrFunc[i] = corrStats[i].GetMean()/var/double(nSamples-i);
}
return corrFunc;
}

double Est_AutoCorr::GetTau() {
// if(nSamples<maxHist) {
// cout << "DES Warning: number of samples less than max history in Est_AutoCorr" << endl;
// cout << nSamples << " < " << maxHist << endl;
// }

// Compute corrFunction
arma::vec logCorr(GetCorr());
logCorr.print("logCorr");
cout << "nSamples = " << nSamples << endl;
//find(logCorr < 0.1, 1).print("first < 0.1")
arma::uvec firstIncrease = find(logCorr < 0.1, 1);
int stopIndex;
if(firstIncrease.size() == 1) {
stopIndex = arma::as_scalar(find(logCorr < 0.1, 1));
}
else {
cout << "Warning: AutoCorr never increases in Est_AutoCorr!" << endl;
stopIndex = min(nSamples,maxHist);
}
cout << "Below 0.1 at " << stopIndex << endl;
bool found=false;
int i = 0;
while(!found && i<corrStats.size()-1) {
i++;
found = logCorr(i-1) < logCorr(i);
}
cout << "first turn up at " << i << endl;
if(i < stopIndex) {
if(logCorr(i) < 0.3) {
stopIndex = i;
}
else {
cout << "Warning: Autocorrelation function turns up while above 0.3 in Est_AutoCorr.cpp" << endl;
stopIndex = (stopIndex+i)/2;
}
}

//stopIndex += i;
//stopIndex /= 2;
cout << "Final stop index = " << stopIndex << endl;
logCorr = logCorr.subvec(0,stopIndex);
int nHist = logCorr.size();
logCorr = log(abs(logCorr));
arma::vec tStep = logCorr;
for(int i=0; i<nHist; i++) {
tStep(i) = double(i);
}
tStep.print("index");
logCorr.print("logCorr");
double tau = -1.0/arma::as_scalar(tStep.t()*logCorr)*double(nHist*(nHist-1)*(2*nHist-1))/6.0;
cout << "tau = " << tau << endl;

exit(-1);

return tau;
}
36 changes: 36 additions & 0 deletions Est_AutoCorr.h
Original file line number Diff line number Diff line change
@@ -0,0 +1,36 @@
/*
* Est_AutoCorr.h
*
* Created on: March 3, 2015
* Author: dstuck
*/

#ifndef EST_AUTOCORR_H_
#define EST_AUTOCORR_H_

#include<cmath>
#include<deque>
#include<vector>
#include "armadillo"
#include "Stats.h"

using namespace std;

class Est_AutoCorr {
public:
Est_AutoCorr();
Est_AutoCorr(int);
virtual ~Est_AutoCorr();
void Reset();
void AddVal(double);
vector<double> GetCorr();
double GetTau();

int maxHist;
int nSamples;
deque<double> hist;
Stats totalStats;
vector<Stats> corrStats;
};

#endif /* EST_AUTOCORR_H_ */
14 changes: 9 additions & 5 deletions Makefile
Original file line number Diff line number Diff line change
@@ -1,13 +1,17 @@
#CXX=icpc
#CXX=g++
#CXXFLAGS=-O0 -g $(INCLUDES)
CXXFLAGS=-O0 -g -DDEBUG
F77FLAGS=-cxxlib -nofor-main
CSOURCES := $(wildcard *.cpp)
ARMAFLAGS=-L/Users/dstuck/tools/armadillo -I/Users/dstuck/tools/armadillo/include -larmadillo -framework Accelerate
#CXXFLAGS=-O0 -g -DDEBUG
#F77FLAGS=-cxxlib -nofor-main
CXXFLAGS=-O0 -g -DDEBUG $(ARMAFLAGS)
F77FLAGS=-cxxlib -nofor-main -Xlinker -no_compact_unwind -framework Accelerate
#CSOURCES := $(wildcard *.cpp)
CSOURCES := $(wildcard *.cpp,optking/*.cpp)
COBJECTS := $(patsubst %.cpp,%.o,$(CSOURCES))
CDEPFILES := $(patsubst %.cpp,%.d,$(CSOURCES))
#LIBTINKER := -L tinker/libtinker.a
LIBTINKER := tinker/active.o tinker/analysis.o tinker/angles.o tinker/attach.o tinker/basefile.o tinker/beeman.o tinker/bicubic.o tinker/bitors.o tinker/bonds.o tinker/born.o tinker/bounds.o tinker/bussi.o tinker/calendar.o tinker/center.o tinker/chkpole.o tinker/chkring.o tinker/chkxyz.o tinker/cholesky.o tinker/clock.o tinker/cluster.o tinker/column.o tinker/command.o tinker/connect.o tinker/connolly.o tinker/control.o tinker/cspline.o tinker/cutoffs.o tinker/deflate.o tinker/delete.o tinker/diagq.o tinker/diffeq.o tinker/eangang.o tinker/eangang1.o tinker/eangang2.o tinker/eangang3.o tinker/eangle.o tinker/eangle1.o tinker/eangle2.o tinker/eangle3.o tinker/ebond.o tinker/ebond1.o tinker/ebond2.o tinker/ebond3.o tinker/ebuck.o tinker/ebuck1.o tinker/ebuck2.o tinker/ebuck3.o tinker/echarge.o tinker/echarge1.o tinker/echarge2.o tinker/echarge3.o tinker/echgdpl.o tinker/echgdpl1.o tinker/echgdpl2.o tinker/echgdpl3.o tinker/edipole.o tinker/edipole1.o tinker/edipole2.o tinker/edipole3.o tinker/egauss.o tinker/egauss1.o tinker/egauss2.o tinker/egauss3.o tinker/egeom.o tinker/egeom1.o tinker/egeom2.o tinker/egeom3.o tinker/ehal.o tinker/ehal1.o tinker/ehal2.o tinker/ehal3.o tinker/eimprop.o tinker/eimprop1.o tinker/eimprop2.o tinker/eimprop3.o tinker/eimptor.o tinker/eimptor1.o tinker/eimptor2.o tinker/eimptor3.o tinker/elj.o tinker/elj1.o tinker/elj2.o tinker/elj3.o tinker/embed.o tinker/emetal.o tinker/emetal1.o tinker/emetal2.o tinker/emetal3.o tinker/emm3hb.o tinker/emm3hb1.o tinker/emm3hb2.o tinker/emm3hb3.o tinker/empole.o tinker/empole1.o tinker/empole2.o tinker/empole3.o tinker/energy.o tinker/dstuckVibrate.o tinker/dstuckEnergy.o tinker/eopbend.o tinker/eopbend1.o tinker/eopbend2.o tinker/eopbend3.o tinker/eopdist.o tinker/eopdist1.o tinker/eopdist2.o tinker/eopdist3.o tinker/epitors.o tinker/epitors1.o tinker/epitors2.o tinker/epitors3.o tinker/erf.o tinker/erxnfld.o tinker/erxnfld1.o tinker/erxnfld2.o tinker/erxnfld3.o tinker/esolv.o tinker/esolv1.o tinker/esolv2.o tinker/esolv3.o tinker/estrbnd.o tinker/estrbnd1.o tinker/estrbnd2.o tinker/estrbnd3.o tinker/estrtor.o tinker/estrtor1.o tinker/estrtor2.o tinker/estrtor3.o tinker/etors.o tinker/etors1.o tinker/etors2.o tinker/etors3.o tinker/etortor.o tinker/etortor1.o tinker/etortor2.o tinker/etortor3.o tinker/eurey.o tinker/eurey1.o tinker/eurey2.o tinker/eurey3.o tinker/evcorr.o tinker/extra.o tinker/extra1.o tinker/extra2.o tinker/extra3.o tinker/fatal.o tinker/fft3d.o tinker/fftpack.o tinker/field.o tinker/final.o tinker/flatten.o tinker/freeunit.o tinker/geometry.o tinker/getint.o tinker/getkey.o tinker/getmol.o tinker/getmol2.o tinker/getnumb.o tinker/getpdb.o tinker/getprm.o tinker/getref.o tinker/getstring.o tinker/gettext.o tinker/getword.o tinker/getxyz.o tinker/ghmcstep.o tinker/gradient.o tinker/gradrgd.o tinker/gradrot.o tinker/groups.o tinker/grpline.o tinker/gyrate.o tinker/hessian.o tinker/hessrgd.o tinker/hessrot.o tinker/hybrid.o tinker/image.o tinker/impose.o tinker/induce.o tinker/inertia.o tinker/initatom.o tinker/initial.o tinker/initprm.o tinker/initres.o tinker/initrot.o tinker/insert.o tinker/invbeta.o tinker/invert.o tinker/jacobi.o tinker/kangang.o tinker/kangle.o tinker/katom.o tinker/kbond.o tinker/kcharge.o tinker/kdipole.o tinker/kewald.o tinker/kgeom.o tinker/kimprop.o tinker/kimptor.o tinker/kinetic.o tinker/kmetal.o tinker/kmpole.o tinker/kopbend.o tinker/kopdist.o tinker/korbit.o tinker/kpitors.o tinker/kpolar.o tinker/ksolv.o tinker/kstrbnd.o tinker/kstrtor.o tinker/ktors.o tinker/ktortor.o tinker/kurey.o tinker/kvdw.o tinker/lattice.o tinker/lbfgs.o tinker/lights.o tinker/makeint.o tinker/makeref.o tinker/makexyz.o tinker/maxwell.o tinker/mdinit.o tinker/mdrest.o tinker/mdsave.o tinker/mdstat.o tinker/mechanic.o tinker/merge.o tinker/molecule.o tinker/moments.o tinker/mutate.o tinker/nblist.o tinker/nextarg.o tinker/nexttext.o tinker/nose.o tinker/nspline.o tinker/number.o tinker/numeral.o tinker/numgrad.o tinker/ocvm.o tinker/openend.o tinker/optsave.o tinker/orbital.o tinker/orient.o tinker/orthog.o tinker/overlap.o tinker/picalc.o tinker/pmestuff.o tinker/pmpb.o tinker/polymer.o tinker/precise.o tinker/pressure.o tinker/prmkey.o tinker/promo.o tinker/prtdyn.o tinker/prterr.o tinker/prtint.o tinker/prtmol2.o tinker/prtpdb.o tinker/prtprm.o tinker/prtseq.o tinker/prtxyz.o tinker/quatfit.o tinker/random.o tinker/rattle.o tinker/readdyn.o tinker/readgau.o tinker/readint.o tinker/readmol.o tinker/readmol2.o tinker/readpdb.o tinker/readprm.o tinker/readseq.o tinker/readxyz.o tinker/replica.o tinker/respa.o tinker/rgdstep.o tinker/rings.o tinker/rmsfit.o tinker/rotlist.o tinker/rotpole.o tinker/sdstep.o tinker/search.o tinker/server.o tinker/shakeup.o tinker/sigmoid.o tinker/sktstuff.o tinker/sort.o tinker/square.o tinker/suffix.o tinker/surface.o tinker/surfatom.o tinker/switch.o tinker/temper.o tinker/tncg.o tinker/torphase.o tinker/torque.o tinker/torsions.o tinker/trimtext.o tinker/unitcell.o tinker/verlet.o tinker/version.o tinker/volume.o tinker/xyzatm.o tinker/zatom.o
LIBTINKER := tinker/active.o tinker/analysis.o tinker/angles.o tinker/attach.o tinker/basefile.o tinker/beeman.o tinker/bicubic.o tinker/bitors.o tinker/bonds.o tinker/born.o tinker/bounds.o tinker/bussi.o tinker/calendar.o tinker/center.o tinker/chkpole.o tinker/chkring.o tinker/chkxyz.o tinker/cholesky.o tinker/clock.o tinker/cluster.o tinker/column.o tinker/command.o tinker/connect.o tinker/connolly.o tinker/control.o tinker/cspline.o tinker/cutoffs.o tinker/deflate.o tinker/delete.o tinker/diagq.o tinker/diffeq.o tinker/eangang.o tinker/eangang1.o tinker/eangang2.o tinker/eangang3.o tinker/eangle.o tinker/eangle1.o tinker/eangle2.o tinker/eangle3.o tinker/ebond.o tinker/ebond1.o tinker/ebond2.o tinker/ebond3.o tinker/ebuck.o tinker/ebuck1.o tinker/ebuck2.o tinker/ebuck3.o tinker/echarge.o tinker/echarge1.o tinker/echarge2.o tinker/echarge3.o tinker/echgdpl.o tinker/echgdpl1.o tinker/echgdpl2.o tinker/echgdpl3.o tinker/edipole.o tinker/edipole1.o tinker/edipole2.o tinker/edipole3.o tinker/egauss.o tinker/egauss1.o tinker/egauss2.o tinker/egauss3.o tinker/egeom.o tinker/egeom1.o tinker/egeom2.o tinker/egeom3.o tinker/ehal.o tinker/ehal1.o tinker/ehal2.o tinker/ehal3.o tinker/eimprop.o tinker/eimprop1.o tinker/eimprop2.o tinker/eimprop3.o tinker/eimptor.o tinker/eimptor1.o tinker/eimptor2.o tinker/eimptor3.o tinker/elj.o tinker/elj1.o tinker/elj2.o tinker/elj3.o tinker/embed.o tinker/emetal.o tinker/emetal1.o tinker/emetal2.o tinker/emetal3.o tinker/emm3hb.o tinker/emm3hb1.o tinker/emm3hb2.o tinker/emm3hb3.o tinker/empole.o tinker/empole1.o tinker/empole2.o tinker/empole3.o tinker/energy.o tinker/dstuckVibrate.o tinker/dstuckOptimize.o tinker/dstuckEnergy.o tinker/eopbend.o tinker/eopbend1.o tinker/eopbend2.o tinker/eopbend3.o tinker/eopdist.o tinker/eopdist1.o tinker/eopdist2.o tinker/eopdist3.o tinker/epitors.o tinker/epitors1.o tinker/epitors2.o tinker/epitors3.o tinker/erf.o tinker/erxnfld.o tinker/erxnfld1.o tinker/erxnfld2.o tinker/erxnfld3.o tinker/esolv.o tinker/esolv1.o tinker/esolv2.o tinker/esolv3.o tinker/estrbnd.o tinker/estrbnd1.o tinker/estrbnd2.o tinker/estrbnd3.o tinker/estrtor.o tinker/estrtor1.o tinker/estrtor2.o tinker/estrtor3.o tinker/etors.o tinker/etors1.o tinker/etors2.o tinker/etors3.o tinker/etortor.o tinker/etortor1.o tinker/etortor2.o tinker/etortor3.o tinker/eurey.o tinker/eurey1.o tinker/eurey2.o tinker/eurey3.o tinker/evcorr.o tinker/extra.o tinker/extra1.o tinker/extra2.o tinker/extra3.o tinker/fatal.o tinker/fft3d.o tinker/fftpack.o tinker/field.o tinker/final.o tinker/flatten.o tinker/freeunit.o tinker/geometry.o tinker/getint.o tinker/getkey.o tinker/getmol.o tinker/getmol2.o tinker/getnumb.o tinker/getpdb.o tinker/getprm.o tinker/getref.o tinker/getstring.o tinker/gettext.o tinker/getword.o tinker/getxyz.o tinker/ghmcstep.o tinker/gradient.o tinker/gradrgd.o tinker/gradrot.o tinker/groups.o tinker/grpline.o tinker/gyrate.o tinker/hessian.o tinker/hessrgd.o tinker/hessrot.o tinker/hybrid.o tinker/image.o tinker/impose.o tinker/induce.o tinker/inertia.o tinker/initatom.o tinker/initial.o tinker/initprm.o tinker/initres.o tinker/initrot.o tinker/insert.o tinker/invbeta.o tinker/invert.o tinker/jacobi.o tinker/kangang.o tinker/kangle.o tinker/katom.o tinker/kbond.o tinker/kcharge.o tinker/kdipole.o tinker/kewald.o tinker/kgeom.o tinker/kimprop.o tinker/kimptor.o tinker/kinetic.o tinker/kmetal.o tinker/kmpole.o tinker/kopbend.o tinker/kopdist.o tinker/korbit.o tinker/kpitors.o tinker/kpolar.o tinker/ksolv.o tinker/kstrbnd.o tinker/kstrtor.o tinker/ktors.o tinker/ktortor.o tinker/kurey.o tinker/kvdw.o tinker/lattice.o tinker/lbfgs.o tinker/lights.o tinker/makeint.o tinker/makeref.o tinker/makexyz.o tinker/maxwell.o tinker/mdinit.o tinker/mdrest.o tinker/mdsave.o tinker/mdstat.o tinker/mechanic.o tinker/merge.o tinker/molecule.o tinker/moments.o tinker/mutate.o tinker/nblist.o tinker/nextarg.o tinker/nexttext.o tinker/nose.o tinker/nspline.o tinker/number.o tinker/numeral.o tinker/numgrad.o tinker/ocvm.o tinker/openend.o tinker/optsave.o tinker/orbital.o tinker/orient.o tinker/orthog.o tinker/overlap.o tinker/picalc.o tinker/pmestuff.o tinker/pmpb.o tinker/polymer.o tinker/precise.o tinker/pressure.o tinker/prmkey.o tinker/promo.o tinker/prtdyn.o tinker/prterr.o tinker/prtint.o tinker/prtmol2.o tinker/prtpdb.o tinker/prtprm.o tinker/prtseq.o tinker/prtxyz.o tinker/quatfit.o tinker/random.o tinker/rattle.o tinker/readdyn.o tinker/readgau.o tinker/readint.o tinker/readmol.o tinker/readmol2.o tinker/readpdb.o tinker/readprm.o tinker/readseq.o tinker/readxyz.o tinker/replica.o tinker/respa.o tinker/rgdstep.o tinker/rings.o tinker/rmsfit.o tinker/rotlist.o tinker/rotpole.o tinker/sdstep.o tinker/search.o tinker/server.o tinker/shakeup.o tinker/sigmoid.o tinker/sktstuff.o tinker/sort.o tinker/square.o tinker/suffix.o tinker/surface.o tinker/surfatom.o tinker/switch.o tinker/temper.o tinker/tncg.o tinker/torphase.o tinker/torque.o tinker/torsions.o tinker/trimtext.o tinker/unitcell.o tinker/verlet.o tinker/version.o tinker/volume.o tinker/xyzatm.o tinker/zatom.o

# LDFLAGS=-L/opt/intel/Compiler/11.1/080/Frameworks/mkl/lib/em64t/ -lblas -L/usr/lib/libshell/ -lshell
.PHONY: clean
Expand All @@ -27,6 +31,6 @@ clean:

$(CDEPFILES): %.d: %.cpp
set -e; rm -f $@; \
$(CXX) -M $< > $@.$$$$; \
$(CXX) -M $(ARMAFLAGS) $< > $@.$$$$; \
sed 's,\($*\)\.o[ :]*,\1.o $@ : ,g' < $@.$$$$ > $@; \
rm -f $@.$$$$
3 changes: 2 additions & 1 deletion PhysicsUtil.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -10,7 +10,8 @@
PhysicsUtil::PhysicsUtil() {
// Makes sure it gets set
numInit = -1;
numFrozModes = 0;
lowFrozModes = 0;
highFrozModes = 0;
lambdaTI = 1.0;
deltaAbInit = false;
charge = 0;
Expand Down
3 changes: 2 additions & 1 deletion PhysicsUtil.h
Original file line number Diff line number Diff line change
Expand Up @@ -28,7 +28,8 @@ class PhysicsUtil {
double morseDE;
double morseAlpha;
int numInit;
int numFrozModes;
int lowFrozModes;
int highFrozModes;
double lambdaTI;
bool deltaAbInit;
int charge;
Expand Down
11 changes: 8 additions & 3 deletions Rho_HO.cpp
Original file line number Diff line number Diff line change
Expand Up @@ -7,7 +7,7 @@

#include "Rho_HO.h"

Rho_HO::Rho_HO(vector<double> w, int nFrozModes) : omega(w), numFrozModes(nFrozModes) {
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) {
Expand Down Expand Up @@ -62,12 +62,17 @@ double Rho_HO::Estimate(vector<Particle> slice1, vector<Particle> slice2, double
}

double est = 0;
for(int j=0; j< numFrozModes; j++) {
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=numFrozModes; j<(int)slice1.size(); 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]);
Expand Down
5 changes: 3 additions & 2 deletions Rho_HO.h
Original file line number Diff line number Diff line change
Expand Up @@ -16,7 +16,7 @@ using namespace std;

class Rho_HO: public Propagator {
public:
Rho_HO(vector<double>, int);
Rho_HO(vector<double>, int, int=0);
Rho_HO(double, int);
virtual ~Rho_HO();

Expand All @@ -29,7 +29,8 @@ class Rho_HO: public Propagator {
string GetType();

vector<double> omega;
int numFrozModes;
int lowFrozModes;
int highFrozModes;
};

#endif /* RHOHO_H_ */
Loading

0 comments on commit 0880537

Please sign in to comment.