main.cpp 6.58 KB
Newer Older
Florian Goth's avatar
Florian Goth committed
1
2
#include <array>
#include <vector>
3
#include <iostream>
4
#include <string>
5
#include <cstdlib>
6
#include <fstream>
Florian Goth's avatar
Florian Goth committed
7
#include "rndwrapper.h"
8
#include "regular_lattice.h"
Manuel Schrauth's avatar
Manuel Schrauth committed
9
#include "vectorhelpers.h"
Florian Goth's avatar
Florian Goth committed
10
#include "registry.h"
Florian Goth's avatar
Florian Goth committed
11

12
13
14
15
using std::cout;
using std::endl;
using std::flush;
using std::ofstream;
Florian Goth's avatar
Florian Goth committed
16

17
#include "systemtools.h"
Florian Goth's avatar
Florian Goth committed
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
//two examples on how to extend the parsing capapbilities of the registry
template <typename T>
struct GetTrait<std::vector<T> >//helper trait to break up a string at various predefined seperators
{
    static std::vector<T> Convert(std::string& arg)
    {
        const std::string delim("; ,");
        std::vector<T> retval;
        std::size_t pos = 0;
        std::size_t posold = 0;
        do
        {
            retval.push_back(GetTrait<T>::Convert(arg.substr(posold, ( pos = arg.find_first_of(delim, posold) ) - posold) ));
        } while ((posold = arg.find_first_not_of(delim, pos) ) != std::string::npos);
        return retval;
    }
};

template <>
struct GetTrait<std::vector<std::string> >//helper trait to break up a string at various predefined seperators
{
    static std::vector<std::string> Convert(std::string arg)
    {
        const std::string delim("; ,");
        std::vector<std::string> retval;
        std::size_t pos = 0;
        std::size_t posold = 0;
        do
        {
            retval.push_back(arg.substr(posold, ( pos = arg.find_first_of(delim, posold) ) - posold) );
        } while ((posold = arg.find_first_not_of(delim, pos) ) != std::string::npos);
        return retval;
    }
};
Florian Goth's avatar
Florian Goth committed
53

54
55
class Grid 
{
Florian Goth's avatar
Florian Goth committed
56
57
    std::vector<int> getnbr(int i);
};
Florian Goth's avatar
Florian Goth committed
58

59
60


Florian Goth's avatar
Florian Goth committed
61
62
template <class StateVector>
class Interaction
Florian Goth's avatar
Florian Goth committed
63
{
64
65
	public:
	    double J;
66
	    virtual StateVector operator() (const StateVector& phi_i) = 0;
Florian Goth's avatar
Florian Goth committed
67
68
};

69

Florian Goth's avatar
Florian Goth committed
70
template <class StateVector, typename CouplingType>
Florian Goth's avatar
Florian Goth committed
71
class OnSite
Florian Goth's avatar
Florian Goth committed
72
{
73
	public: 
Florian Goth's avatar
Florian Goth committed
74
	    CouplingType h;
75
	    virtual CouplingType operator() (const StateVector& phi) = 0;
76
	private:
Florian Goth's avatar
Florian Goth committed
77
78
};

79
80
81


template <class StateSpace, class StateVector>
Florian Goth's avatar
Florian Goth committed
82
class MultiSite
Florian Goth's avatar
Florian Goth committed
83
{
84
85
	public: 
	    const double k;
86
	    virtual double operator() (const StateVector& sv, int svpos, StateSpace s) = 0;
87
	private:
Florian Goth's avatar
Florian Goth committed
88
89
};

90
91
92
93
94


// ------- elementary state vector calculus

template <class StateVector>
95
StateVector operator + (StateVector lhs,  StateVector rhs)
Florian Goth's avatar
Florian Goth committed
96
{
Florian Goth's avatar
bugfix    
Florian Goth committed
97
98
99
    StateVector res(lhs);
    for(int i = 0; i < std::tuple_size<StateVector>::value; ++i)
    res[i] += rhs[i];
Florian Goth's avatar
Florian Goth committed
100
101
102
    return res;
}

103
template <class StateVector>
104
StateVector operator - (StateVector lhs,  StateVector rhs)
Florian Goth's avatar
Florian Goth committed
105
{
Florian Goth's avatar
bugfix    
Florian Goth committed
106
107
108
    StateVector res(lhs);
    for(int i = 0; i < std::tuple_size<StateVector>::value; ++i)
    res[i] -= rhs[i];
Florian Goth's avatar
Florian Goth committed
109
110
111
    return res;
}

Florian Goth's avatar
Florian Goth committed
112
113
114
115
116
inline double dot(const double& a, const double& b)
{
    return a*b;
}

117
118
template<class VecType>
inline typename VecType::value_type dot(const VecType& a, const VecType& b)
Florian Goth's avatar
Florian Goth committed
119
{
Manuel Schrauth's avatar
Manuel Schrauth committed
120
    typedef typename VecType::value_type FPType;
121
    return std::inner_product(begin(a), end(a), begin(b), 0.0);
Florian Goth's avatar
Florian Goth committed
122
123
}

Florian Goth's avatar
Florian Goth committed
124

125
126
127
128
129
130
131
132
133
134
template <class StateVector>
inline void reflect(StateVector& vec, const StateVector mirror)
{
	const int SymD = std::tuple_size<StateVector>::value;
	
	const double dotp = dot(vec,mirror);

	for (int i=0; i<SymD; i++) vec[i] -= 2*dotp*mirror[i];
}	

135
136
137
138
139
140
141
template <class Container>
inline void normalize(Container& a)
{
	typename Container::value_type tmp_abs=std::sqrt(dot(a, a));

	for (int i = 0; i < a.size(); ++i) a[i] /= tmp_abs;
}
142
143
144



145
146
147
148
149
150
151
152
153
154
template <class StateVector>
inline void coutsv(StateVector& vec)
{
	const int SymD = std::tuple_size<StateVector>::value;
	
	for (int i=0; i<SymD; i++) cout << vec[i] << "\t";
	cout << endl;
}	


155
// ---------------------------------------
Florian Goth's avatar
Florian Goth committed
156

157
158
const int myid = 0; // remove once a parallelization is available

159
160
#include "Heisenberg.h"
#include "Ising.h"
Manuel Schrauth's avatar
Manuel Schrauth committed
161
#include "Phi4.h"
162
#include "BlumeCapel.h"
163
#include "XXZAntiferro.h"
164

165
#include "marqov.h"
166

Florian Goth's avatar
Florian Goth committed
167
168
int main()
{
169
	// read config files
Manuel Schrauth's avatar
Manuel Schrauth committed
170
	RegistryDB registry("../src/config");
171

172

Manuel Schrauth's avatar
Manuel Schrauth committed
173
	// remove old output and prepare new one
174
175
176
177
178
	std::string outdir = registry.Get<std::string>("mc", "IO", "outdir" );
	std::string logdir = registry.Get<std::string>("mc", "IO", "logdir" );

	std::string command;
	command = "rm -r " + outdir;
Manuel Schrauth's avatar
Manuel Schrauth committed
179
	system(command.c_str());
180
181
182
	command = "rm -r " + logdir;
	system(command.c_str());

Manuel Schrauth's avatar
Manuel Schrauth committed
183
	makeDir(outdir);
184
	makeDir(logdir);
185
186
	

Manuel Schrauth's avatar
Manuel Schrauth committed
187
188
	// ------------------ live view -----------------------
	/*
Manuel Schrauth's avatar
Manuel Schrauth committed
189

190
191
192
		const int L = 30;
		RegularLattice lattice(L, 3);
		std::string outfile = outdir+"temp.h5";
193
		Marqov<RegularLattice, XXZAntiferro<double,double> > marqov(lattice, outfile, 1/0.36);
194
		marqov.init_hot();
195
		const int ncluster = 0;
196
197
198
		const int nsweeps  = L/2; 
		marqov.wrmploop(50, ncluster, nsweeps);
		marqov.gameloop_liveview();
199

200
	*/
Manuel Schrauth's avatar
Manuel Schrauth committed
201
	// ----------------------------------------------------
Manuel Schrauth's avatar
Manuel Schrauth committed
202

203

204
205
	// extract Monte Carlo parameters
	auto  nL = registry.Get<std::vector<int> >("mc", "General", "nL" );
206
207
208
	int    nbeta     = registry.Get<int>("mc", "General", "nbeta" );
	double betastart = registry.Get<double>("mc", "General", "betastart" );
	double betaend   = registry.Get<double>("mc", "General", "betaend" );
Manuel Schrauth's avatar
Manuel Schrauth committed
209
210
	double betastep = (betaend - betastart) / double(nbeta);

211

212
213
214
	// write temperatures in logfile
	std::string logfile = registry.Get<std::string>("mc", "IO", "logfile" );
	std::ofstream os(logdir+"/"+logfile);
Manuel Schrauth's avatar
Manuel Schrauth committed
215
216
217
	for (int i=0; i<nbeta; i++) os << betastart + i*betastep << endl;
	os.close();

Manuel Schrauth's avatar
Manuel Schrauth committed
218
219
220
221
222
223


	// lattice dimension
	const int dim = 3;

	// lattice size loop
224
	for (int j=0; j<nL.size(); j++)
Manuel Schrauth's avatar
Manuel Schrauth committed
225
	{
Manuel Schrauth's avatar
Manuel Schrauth committed
226
		// extract lattice size and prepare directories
Manuel Schrauth's avatar
Manuel Schrauth committed
227
228
		const int L = nL[j];
		cout << endl << "L = " << L << endl << endl;
229
		cout << outdir+std::to_string(L) << endl;
Manuel Schrauth's avatar
Manuel Schrauth committed
230
		makeDir(outdir+std::to_string(L));
231

Manuel Schrauth's avatar
Manuel Schrauth committed
232
		// temperature loop
233
234
235
236
237
		for (int i=0; i<nbeta; i++)
		{
			double currentbeta = betastart + i*betastep; 
			cout << "beta = " << currentbeta << endl;

Manuel Schrauth's avatar
Manuel Schrauth committed
238
239
			// set up lattice
			RegularLattice lattice(L, dim);
240
		
Manuel Schrauth's avatar
Manuel Schrauth committed
241
			// set up outfile
242
			std::string outfile = outdir+"/"+std::to_string(L)+"/"+std::to_string(i)+".h5";
Manuel Schrauth's avatar
Manuel Schrauth committed
243
244
245
246
247
248

			// set up model
//			Marqov<RegularLattice, Ising<int> > marqov(lattice, currentbeta, outfile);
//			Marqov<RegularLattice, BlumeCapel<int> > marqov(lattice, currentbeta, outfile);
//			Marqov<RegularLattice, Heisenberg<double,double> > marqov(lattice, outfile, registry);
//			Marqov<RegularLattice, Phi4<double,double> > marqov(lattice, outfile, currentbeta);
249
			Marqov<RegularLattice, XXZAntiferro<double,double> > marqov(lattice, outfile, currentbeta);
250

251

Manuel Schrauth's avatar
Manuel Schrauth committed
252
			// number of cluster updates and metropolis sweeps
253
254
			const int ncluster = 0;
			const int nsweeps  = L; 
Manuel Schrauth's avatar
Manuel Schrauth committed
255

256

Manuel Schrauth's avatar
Manuel Schrauth committed
257
			// number of EMCS during relaxation and measurement
258
259
			const int nrlx = 2500;
			const int nmsr = 5000;  
Manuel Schrauth's avatar
Manuel Schrauth committed
260
261


Manuel Schrauth's avatar
Manuel Schrauth committed
262
263
264
			// perform simulation
			marqov.init_hot();
			marqov.wrmploop(nrlx, ncluster, nsweeps);
Manuel Schrauth's avatar
Manuel Schrauth committed
265
			marqov.gameloop(nmsr, ncluster, nsweeps);
Manuel Schrauth's avatar
Manuel Schrauth committed
266
267
268

		}
	}
Florian Goth's avatar
Florian Goth committed
269
}