Ssh.h 13.2 KB
Newer Older
1
2
3
4
5
6
7
#ifndef SSH_H
#define SSH_H
#include <array>
#include <tuple>
#include <string>
#include <complex>
#include <functional>
8
#include <math.h> // for debug, provides "isnan"
9
10
#include "../hamparts.h"
#include "../metropolis.h"
11

12
13
#define SSH_2D

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

// ------------------------------ OBSERVABLES ---------------------------

class SSHMag
{
	public:
		std::string name, desc;
		template <class StateSpace, class Grid> 
		double measure(const StateSpace& statespace, const Grid& grid) 
		{
			const int N = grid.size();
			double mag = 0.0;

			for (int i=0; i<N; i++) 
			{
				 mag += statespace[i][0];
				}
			return mag/double(N);
		}
		SSHMag() : name("m"), desc("The Magnetization of the SSH Modell") {}
};


class SSHMagSq
{
	public:
		std::string name, desc;
		template <class StateSpace, class Grid> 
		double measure(const StateSpace& statespace, const Grid& grid) 
		{
			const int N = grid.size();
			double mag = 0.0;

			for (int i=0; i<N; i++) 
			{
				 mag += pow(statespace[i][0],2);
				}
			return mag/double(N);
		}

		SSHMagSq() : name("msq"), desc("...") {}
};




class SSHTwoPointCorrSpace
{
	public:
	std::string name, desc;
	template <class StateSpace, class Grid>
65
	std::array<double, 64> measure(const StateSpace& statespace, const Grid& grid)
66
67
68
69
70
71
72
73
74
75
76
//	std::vector<double> measure(const StateSpace& statespace, const Grid& grid) // not yet supported
	{

		const int Ls = grid.len;
		const int Lt = grid.lentime;

		const double norml = 1.0 / Ls / Ls;

//		std::vector<double> retval;
//		for (int i=0; i<10; i++) retval.push_back(0);

77
78
		std::array<double,64> retval;
		for (int i=0; i<64; i++) retval[i] = 0;
79
80
81
82
83

		for (int i=0; i<Lt; i++)
		{
			for (int j=0; j<Ls; j++)
			{
Manuel Schrauth's avatar
Manuel Schrauth committed
84
85
				const int idx1 = i*Ls+j;

86
87
				for (int k=0; k<Ls; k++)
				{
88
//					if (k==j) continue;
89
90
91
92
93

					const int idx2 = i*Ls+k;

					const double x1 = statespace[idx1][0];
					const double x2 = statespace[idx2][0];
Manuel Schrauth's avatar
Manuel Schrauth committed
94

95
96
97
98
99
100
					//// actual coordinates not necessary for 1D
					// const auto c1 = grid.getcrds(idx1)[0];
					// const auto c2 = grid.getcrds(idx2)[0];
					// auto dist = c1-c2;
					// if (fabs(dist) > 0.5) dist = 1.0 - fabs(dist);

Manuel Schrauth's avatar
Manuel Schrauth committed
101
102
					int diff = std::abs(j-k);
					if (diff>Ls/2) diff -= Ls/2;
103
104
105
106
107
108
109

					retval[diff] += norml*x1*x2;
				}
			}
		}

		return retval;
Manuel Schrauth's avatar
Manuel Schrauth committed
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
	}
					
	
	SSHTwoPointCorrSpace() : name("corr"), desc("Spatial two-point correlator") {}

};

template <class StateVector>
class SSH_interaction : public Interaction<StateVector> 
{
public:
	SSH_interaction(double m, double dtau)
	{
		this->J = -m/dtau;
	}
	StateVector get (const StateVector& phi) {return phi;};
};




template <class StateVector>
class SSH_onsite : public OnSite<StateVector, double> 
{
	public:
		SSH_onsite(double m, double k, double dtau)
		{
			this->h = m/dtau + k*dtau/2;
		}
		double get (const StateVector& phi) {return dot(phi,phi);}; 
};



Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
145
146
// =========================== Multi-Site Class =================================

147
148
149
150
template <class StateSpace, class StateVector>
class SSH_multisite
{
	public:
151
		double k, beta, dtau, g;
Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
		int L;
		double* gdat;
		int ntau;
		SSH_multisite(double g, double b, double d, int myL) : k(-0.5*g*g), beta(b), dtau(d), L(myL), g(g)
		{
			ntau = std::round(beta/dtau); // ntau should be nothing else than Ltime
			gdat = new double[ntau*L];

			for(int j = 0; j < L; ++j)
			{
				double k = (j*2)*M_PI/double(L);
				// 1D disperson relation
				const double mu = 0;
				double eps = -2.0*std::cos(k)-mu;
				for(int dt = 0; dt < ntau; ++dt)
				{
Manuel Schrauth's avatar
Manuel Schrauth committed
168
					gdat[dt * L + j] = 0.5*std::exp((-beta/2 + dt*dtau)*eps)/std::cosh(beta*eps/2);
Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
169
170
171
172
				}
			}
		}

173
174
175
176
177
178
179
180
181
182
183
184

		// computes the difference already!
		// can be simplified if "suscept" is symmetric

		template <class Lattice>
		double diff (const int rsite,
					const StateVector& svold, 
					const StateVector& svnew, 
					std::vector<int>& nbrs, 
					StateSpace& s,
					Lattice& grid)
		{
185

186
187
188
189
			double retval = 0;
			for (int i=0; i<nbrs.size(); i++)
			{
				auto b1 = svnew-svold;
190
				auto b2 = suscept(grid, rsite, nbrs[i]);
191
192
193
194
195
196
197
198
				auto b3 = s[i];

				auto a1 = dot(b1, mult(b2,b3));
				auto a2 = dot(b3, mult(b2,b1));

				retval = retval + a1 + a2;
			}

199
200
			retval += dot(svnew, mult(suscept(grid, rsite, rsite), svnew));
			retval -= dot(svold, mult(suscept(grid, rsite, rsite), svold));
201
202
203
204

			return retval;
		}

Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
205
		~SSH_multisite() {delete [] gdat;}
Manuel Schrauth's avatar
Manuel Schrauth committed
206

Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
		/**
		The fermi function
		@return An occupation propability
		*/
		template <typename T>
		inline T fermi(const T& e) throw()
		{
		    T xhalf = -e/T(2.0);
		    T retval;
		    if(xhalf > std::log(0.01*std::numeric_limits<T>::max()))
		      retval = 1.0;
		    else
		      if(xhalf < std::log(10*std::numeric_limits<T>::min()))
			retval = 0.0;//FIXME: not yet decided whether setting this to zero is better than setting it to epsilon
			else
			  retval = 0.5*std::exp(xhalf)/std::cosh(xhalf);
		    return retval;
		}
	
226
		
227
		#ifdef SSH_2D
Manuel Schrauth's avatar
Manuel Schrauth committed
228
229
		// Green's function for a square lattice
		// takes two coordinate vectors (x,y,t), that is 2+1 dimensions
Manuel Schrauth's avatar
Manuel Schrauth committed
230
		double green(std::vector<double> c1, std::vector<double> c2)
231
232
233
		{
			std::complex<double> retval = 0;
			std::complex<double> jj(0,1);
Manuel Schrauth's avatar
bufix    
Manuel Schrauth committed
234

Manuel Schrauth's avatar
Manuel Schrauth committed
235
236
237
			// space
			double distx = c1[0]-c2[0];
			double disty = c1[1]-c2[1];
238
239
240
			if (distx < 0) distx = L + distx;
			if (disty < 0) disty = L + disty;

Manuel Schrauth's avatar
Manuel Schrauth committed
241
			// time
Manuel Schrauth's avatar
bufix    
Manuel Schrauth committed
242
243
			const double t1 = c1[2];
			const double t2 = c2[2];
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
			int t1i = floor(t1);
			int t2i = floor(t2);
			int dti = t1i - t2i;
			double signum = 1;
			if (dti < 0) 
			{
				signum = -1; // anti-symmetric behaviour in delta tau
				dti = ntau+dti;
			}

			// compute Greens function in Fourier space
			std::complex<double> dexpkx = std::exp(-2*M_PI/L*jj*distx);
			std::complex<double> dexpky = std::exp(-2*M_PI/L*jj*disty);
			std::complex<double> expkx = 1.0;
			std::complex<double> expky = 1.0;

			for (int jx = 0; jx < L; ++jx)
			{
				for (int jy = 0; jy < L; ++jy)
				{
Manuel Schrauth's avatar
Manuel Schrauth committed
264
					// dispersion relation
265
					double disp = - 2*cos(2*M_PI*jx/L) - 2*cos(2*M_PI*jy/L);
Manuel Schrauth's avatar
Manuel Schrauth committed
266
					// do the summation
267
					retval += expkx * expky * exp(dti*dtau*disp) * fermi(beta*disp);
Manuel Schrauth's avatar
Manuel Schrauth committed
268
269
					// increment Fourier transform
					expky *= dexpky; 
270
271
272
273
				}
				expkx *= dexpkx; // increment Fourier transform
			}

Manuel Schrauth's avatar
Manuel Schrauth committed
274
275
			const double norml = 1.0 / pow(2*L,2); // the number of sites per time slice
			return norml*signum*retval.real();
276
277
		}

278
279
		#else

280

Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
281
		// Green's function for a regular hypercubic lattice in 1D (i.e. a chain)
Manuel Schrauth's avatar
Manuel Schrauth committed
282
		double green(std::vector<double> c1, std::vector<double> c2, int sign1, int sign2)
Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
283
284
285
286
287
		{
			std::complex<double> retval = 0;
			std::complex<double> jj(0,1);
		
			double step = 2*M_PI/double(L);
Manuel Schrauth's avatar
Manuel Schrauth committed
288
		
Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
289
			// spatial coordinates
Manuel Schrauth's avatar
Manuel Schrauth committed
290
291
			const double r1 = c1[0] + sign1*0.5; // account for p.b.c not relevant? (we only use relative distances...)
			const double r2 = c2[0] + sign2*0.5;
Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
292
293
			double dist = r1-r2;
			if (dist < 0) dist = L + dist;
Manuel Schrauth's avatar
Manuel Schrauth committed
294
	
Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
295
296
297
			// temporal coordinates
			const double t1 = c1[1];
			const double t2 = c2[1];
Manuel Schrauth's avatar
Manuel Schrauth committed
298
299
			int t1i = floor(t1);
			int t2i = floor(t2);
Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
300
301
302
303
304
305
306
			int dti = t1i - t2i;
			double signum = 1;
			if (dti < 0) 
			{
				signum = -1; // anti-symmetric behaviour in delta tau
				dti = ntau+dti;
			}
Manuel Schrauth's avatar
Manuel Schrauth committed
307
	
Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
308
			// account for periodic boundaries of the lattice
Manuel Schrauth's avatar
Manuel Schrauth committed
309
			// not needed, says Florian
Manuel Schrauth's avatar
Manuel Schrauth committed
310
311
//			if (dti  > 0.5*ntau) dti = ntau - dti;
//			if (dist > 0.5*L)   dist = L - dist;
Manuel Schrauth's avatar
Manuel Schrauth committed
312
	
Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
313
314
	
			// compute Greens function in Fourier space
Manuel Schrauth's avatar
Manuel Schrauth committed
315
			std::complex<double> dexpk = std::exp(-2*M_PI/L*jj*dist);
Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
316
317
318
319
			std::complex<double> expk = 1.0;
			for (int j = 0; j < L; ++j)
			{
				retval += expk.real()*gdat[dti*L + j];
Manuel Schrauth's avatar
Manuel Schrauth committed
320
				expk   *= dexpk; //loop-carried dependency. breaks vectorization.
Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
321
322
			}
			auto retv = signum*retval.real()/double(L);
Manuel Schrauth's avatar
Manuel Schrauth committed
323

Manuel Schrauth's avatar
tidy up    
Manuel Schrauth committed
324
			return retv;
Manuel Schrauth's avatar
Manuel Schrauth committed
325
		}
326
		#endif
Manuel Schrauth's avatar
Manuel Schrauth committed
327
328
329
330
331
332
333
334
335
336
337
338
339


/*

	 j1/i2      b2     j2/i3      b3      j3/i4      b4
  ------x-----------------x------------------x-----------------x---------- // original model: bonds

  ---------------o-----------------o------------------o----------------    // our model: bonds represented as sites


   e.g.  r(b2(j2)) = r(b3(i3))

*/
340

341
342
343

#ifdef SSH_2D

Manuel Schrauth's avatar
Manuel Schrauth committed
344
	// performs a Wick decomposition
Manuel Schrauth's avatar
Manuel Schrauth committed
345
	std::complex<double> wick(std::vector<double> v0, std::vector<double> v1, std::vector<double> v2, std::vector<double> v3)
Manuel Schrauth's avatar
Manuel Schrauth committed
346
347
	{
		std::complex<double> retval = 0;
348
		if (v0==v1 && v1==v2 && v2==v3 && v3==v0)
Manuel Schrauth's avatar
Manuel Schrauth committed
349
		{
Manuel Schrauth's avatar
Manuel Schrauth committed
350
			retval = green(v0,v0);
Manuel Schrauth's avatar
Manuel Schrauth committed
351
		}
352
		else
Manuel Schrauth's avatar
Manuel Schrauth committed
353
		{
Manuel Schrauth's avatar
Manuel Schrauth committed
354
			retval = green(v0,v1)*green(v2,v3) - green(v0,v3)*green(v2,v1);
Manuel Schrauth's avatar
Manuel Schrauth committed
355

356
		}
Manuel Schrauth's avatar
Manuel Schrauth committed
357
		return retval;
358

Manuel Schrauth's avatar
Manuel Schrauth committed
359
360
	}

361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
#else

  	std::complex<double> wick(std::vector<std::vector<double>> s, std::vector<int> e)
  	{
  		std::complex<double> retval = 0;
  
  		// the square of the number operator of fermions is the number operator
  		if (s[0]==s[1] && s[1]==s[2] && s[2]==s[3] && e[0]==e[1] && e[1]==e[2] && e[2]==e[3])
  		{
  			retval = green(s[0], s[1], e[0], e[1]);
  		}
  
  		// Wick decomposition
  		else 
  		{
  			retval = green(s[0],s[1],e[0],e[1])*green(s[2],s[3],e[2],e[3]) 
  				  - green(s[0],s[3],e[0],e[3])*green(s[2],s[1],e[2],e[1]);
  		}
  
  		return retval;
  	}

#endif


Manuel Schrauth's avatar
Manuel Schrauth committed
386
	template <class Lattice>
387
	double suscept(Lattice& grid, int idx1, int idx2)
Manuel Schrauth's avatar
Manuel Schrauth committed
388
	{
389
		std::complex<double> retval = 0;
390

Manuel Schrauth's avatar
Manuel Schrauth committed
391
392
393
		auto w1 = grid.getcrds(idx1);
		auto w2 = grid.getcrds(idx2);

394
395
396
397
398
399
400
401
402
403
404
405
406
		#ifdef SSH_2D
			const auto t1 = w1[4];
			const auto t2 = w2[4];

//			const auto c1 = wick({w1[0],w1[1],t1}, {w1[2],w1[3],t1}, {w2[0],w2[1],t2}, {w2[2],w2[3],t2});
//			const auto c2 = wick({w1[0],w1[1],t1}, {w1[2],w1[3],t1}, {w2[2],w2[3],t2}, {w2[0],w2[1],t2});
//			const auto c3 = wick({w1[2],w1[3],t1}, {w1[0],w1[1],t1}, {w2[0],w2[1],t2}, {w2[2],w2[3],t2});
//			const auto c4 = wick({w1[2],w1[3],t1}, {w1[0],w1[1],t1}, {w2[2],w2[3],t2}, {w2[0],w2[1],t2});

			const auto c1 = wick({w1[0],w1[2],t1}, {w1[1],w1[3],t1}, {w2[0],w2[2],t2}, {w2[1],w2[3],t2});
			const auto c2 = wick({w1[0],w1[2],t1}, {w1[1],w1[3],t1}, {w2[1],w2[3],t2}, {w2[0],w2[2],t2});
			const auto c3 = wick({w1[1],w1[3],t1}, {w1[0],w1[2],t1}, {w2[0],w2[2],t2}, {w2[1],w2[3],t2});
			const auto c4 = wick({w1[1],w1[3],t1}, {w1[0],w1[2],t1}, {w2[1],w2[3],t2}, {w2[0],w2[2],t2});
407

408
409
//			const std::complex<double> K1 = green({w1[0],w1[1],w1[4]},{w1[2],w1[3],w1[4]}) + green({w1[2],w1[3],w1[4]},{w1[0],w1[1],w1[4]});
//			const std::complex<double> K2 = green({w2[0],w2[1],w2[4]},{w2[2],w2[3],w2[4]}) + green({w2[2],w2[3],w2[4]},{w2[0],w2[1],w2[4]});
410

411
412
413
414
415
416
417
418
419
420
			const std::complex<double> K1 = green({w1[0],w1[2],t1},{w1[1],w1[3],t1}) + green({w1[1],w1[3],t1},{w1[0],w1[2],t1});
			const std::complex<double> K2 = green({w2[0],w2[2],t2},{w2[1],w2[3],t2}) + green({w2[1],w2[3],t2},{w2[0],w2[2],t2});
		#else
			const auto t1 = w1[2];
			const auto t2 = w2[2];

			const auto c1 = wick({w1[0],t1}, {w1[1],t1}, {w2[0],t2}, {w2[1],t2});
			const auto c2 = wick({w1[0],t1}, {w1[1],t1}, {w2[1],t2}, {w2[0],t2});
			const auto c3 = wick({w1[1],t1}, {w1[0],t1}, {w2[0],t2}, {w2[1],t2});
			const auto c4 = wick({w1[1],t1}, {w1[0],t1}, {w2[1],t2}, {w2[0],t2});
421

422
423
424
			const std::complex<double> K1 = green({w1[0],t1},{w1[1],t1}) + green({w1[1],t1},{w1[0],t1});
			const std::complex<double> K2 = green({w2[0],t2},{w2[1],t2}) + green({w2[1],t2},{w2[0],t2});
		#endif
Manuel Schrauth's avatar
Manuel Schrauth committed
425

426
		return std::real(c1+c2+c3+c4-K1*K2);
427

Manuel Schrauth's avatar
Manuel Schrauth committed
428
	}
429

Manuel Schrauth's avatar
Manuel Schrauth committed
430
431
432
433
434
435
436
};






437
438
439
440
441
442
443
444
445
446
447
448
449
450
451

template <class StateVector, class RNG>
class SSH_Initializer
{
	private:
		RNG& rng;
	public:
		SSH_Initializer()   {}
		SSH_Initializer(RNG& rn) : rng(rn) {}

		// specifies how a random new state vector is generated
		// in this case a simple spin flip
		StateVector newsv(const StateVector& svold) 
		{
			StateVector retval(svold); 
Manuel Schrauth's avatar
Manuel Schrauth committed
452
			double amp = 0.7;
453
454
455
456
457
458
459
460
461
462
463
464
465
466
			double diff = rng.real(-amp, amp);
			retval[0] += diff;
			return retval;
		};
};



// ------------------------------ HAMILTONIAN ---------------------------

template <typename SpinType = double>
class SSH
{
	public:
Manuel Schrauth's avatar
Manuel Schrauth committed
467
		double m, k, g, dtau, betaQM, L;
468
469
470
471
472
473
474
		constexpr static int SymD = 1;
		const std::string name;
		typedef std::array<SpinType, SymD> StateVector;
		template <typename RNG>
		using MetroInitializer = SSH_Initializer<StateVector, RNG>;

		static constexpr uint Nalpha = 1;
Manuel Schrauth's avatar
Manuel Schrauth committed
475
		static constexpr uint Nbeta = 1;
Manuel Schrauth's avatar
Manuel Schrauth committed
476
		static constexpr uint Ngamma = 1;
477
		
Manuel Schrauth's avatar
Manuel Schrauth committed
478
		SSH(double m, double k, double g, double bQ, int Ltime, int L) : m(m), k(k), g(g), dtau(bQ/double(Ltime)), L(L), betaQM(bQ), name("SSH")
479
480
		{
			interactions[0] = new SSH_interaction<StateVector>(m, dtau); 
Manuel Schrauth's avatar
Manuel Schrauth committed
481
			onsite[0] = new SSH_onsite<StateVector>(m, k, dtau);
Manuel Schrauth's avatar
Manuel Schrauth committed
482
			multisite[0] = new SSH_multisite<StateVector*,StateVector>(g, betaQM, dtau, L);
483
484
485
486
		}
		
		// instantiate interaction terms (requires pointers)
		Interaction<StateVector>* interactions[Nalpha];
Manuel Schrauth's avatar
Manuel Schrauth committed
487
		OnSite<StateVector, double>* onsite[Nbeta];
Manuel Schrauth's avatar
Manuel Schrauth committed
488
		SSH_multisite<StateVector*,  StateVector>* multisite[Ngamma];
489
490
491
	
		// instantiate and choose observables
		SSHMag       obs_m;
Manuel Schrauth's avatar
Manuel Schrauth committed
492
		SSHMagSq       obs_msq;
Manuel Schrauth's avatar
Manuel Schrauth committed
493
494
		SSHTwoPointCorrSpace obs_corr;
		auto getobs()	{return std::make_tuple(obs_m, obs_msq, obs_corr);}
495
496
497
498
499
500
501
502


		// initialize state space
		template <class StateSpace, class Lattice, class RNG>
		void initstatespace(StateSpace& statespace, Lattice& grid, RNG& rng) const
		{
			for (int i=0; i<grid.size(); i++)
			{
Manuel Schrauth's avatar
Manuel Schrauth committed
503
				if (rng.real() > 0.5) statespace[i][0] = 1;
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
				else statespace[i][0] = -1;
			}
		}


		// using the Wolff cluster algorithm requires to implement
		// the functions 'wolff_coupling' and 'wolff_flip'

		template <class A = bool>
		inline double wolff_coupling(StateVector& sv1, StateVector& sv2, const A a=0) const 
		{
			if (sv1[0] == sv2[0]) return 0.0;
			else return -1.0;
		}

		template <class A = bool>
		inline void wolff_flip(StateVector& sv, const A a=0) const 
		{
			sv[0] *= -1;
		}


};


#endif