To receive notifications about scheduled maintenance, please subscribe to the mailing-list gitlab-operations@sympa.ethz.ch. You can subscribe to the mailing-list at https://sympa.ethz.ch

Commit 21e6e135 authored by Valerio's avatar Valerio
Browse files

Solutions exercise 2

parent a4d8638c
This diff is collapsed.
This diff is collapsed.
all: main data.dat plot
main: trFieldIsingTimeEvolution.cpp
g++ trFieldIsingTimeEvolution.cpp -o main -O3
data.dat: main
./main > data.dat
plot: data.dat
python3 plot.py
import numpy as np
import matplotlib.pyplot as plt
data = np.loadtxt("data.dat")
# plot magnetization
plt.pcolor(data,
vmin=-1.0, vmax=1.0)
plt.xlabel('site')
plt.ylabel('time')
#plt.axis([0,l,tstart,tend])
plt.title("Magnetisation")
plt.colorbar()
plt.show()
#include <valarray>
#include <vector>
#include <iostream>
#include <cmath>
#include <complex>
typedef std::valarray< std::complex<double> > state_vector_t;
typedef std::vector<double> param_vector_t;
typedef unsigned state_t;
const std::complex<double> I(0.,1.);
template <class Vector>
void tstep_transv(unsigned l, const param_vector_t& cosinesPrecomputed, const Vector& sinesPrecomputed, Vector& x, Vector& tmp)
{
state_t dim = 1 << l;
tmp *= 0;
for( int r = 0; r < l ; ++r )
{
// diagonal
for( state_t s = 0; s < dim; ++s )
tmp[s] = cosinesPrecomputed[r] * x[s];
// off-diagonal
for( state_t s = 0; s < dim; ++s )
tmp[s^(1<<r)] += sinesPrecomputed[r] * x[s];
std::swap(x,tmp);
}
}
param_vector_t precomputeCosines(const param_vector_t &hx, double dt){
param_vector_t res(hx.size(),double(0.0));
for(int r = 0; r<hx.size(); ++r){
res[r] = std::cos(dt*hx[r]);
}
return res;
}
template <class Vector>
Vector precomputeSines(const param_vector_t &hx, double dt){
Vector res(std::complex<double>(0.0,0.0), hx.size());
for(int r = 0; r<hx.size(); ++r){
res[r] = I * std::sin(dt*hx[r]);
}
return res;
}
template <class Vector>
void tstep_diag(unsigned l, const param_vector_t &j, double dt, Vector& x)
{
state_t dim = 1 << l;
for( state_t s = 0; s < dim; ++s )
{
double jtotal = 0.;
// +J for parallel, -J for antiparallel neighbors
for( int r = 0; r < l - 1; ++r )
{
jtotal += ((s >> r)^(s >> (r+1)))&1 ? -j[r] : j[r]; // check if spins are parallel at site r and r+1
}
x[s] *= std::exp(I*jtotal*dt);
}
}
template <class Vector>
void evolve(unsigned l, const param_vector_t &j, const param_vector_t &hx, double dt, unsigned n, Vector& x)
{
state_t dim = 1 << l;
state_vector_t tmp(x.size()); // memory for tstep_transv
param_vector_t cosinesPrecomputed = precomputeCosines(hx,dt);
Vector sinesPrecomputed = precomputeSines<Vector>(hx,dt);
tstep_diag(l,j,dt/2.,x);
for(int i=0 ; i<n-1 ; ++i){
tstep_transv(l,cosinesPrecomputed,sinesPrecomputed,x,tmp);
tstep_diag(l,j,dt,x);
}
tstep_transv(l,cosinesPrecomputed,sinesPrecomputed,x,tmp);
tstep_diag(l,j,dt/2.,x);
}
template <class Vector>
void printMagnetization(unsigned l, const Vector& v, std::ostream& outputStream)
{
state_t dim = 1 << l;
for(int r=0; r<l; r++){
double m(0);
for(long s=0; s<dim; s++){
m += std::norm(v[s]) * 2. * ( (bool)(s&(1 << r)) -0.5);
}
outputStream << m << "\t\t";
}
outputStream << std::endl;
}
int main(){
// ======== input
unsigned l = 21; // length of chain
param_vector_t j(l-1,1.0); // Ising coupling
param_vector_t hx(l,0.4); // transverse magnetic field
double tmax = 20.; // maximum time reached in time evolutin
double tmeas = 0.5; // time after which a measurment is performed
double dt = 0.1; // time step used for the evolution
unsigned nmeas = tmax/tmeas; // number of measurements
unsigned nstep = tmeas/dt; // number of steps between measurements
// ======== input
// Starting configuration
state_t initialInd = 1<<10;
state_vector_t x(0., 1 << l);
x[initialInd] = 1;
printMagnetization(l,x,std::cout);
// Do time evolution
for(int n=1; n<=nmeas; ++n){
evolve(l, j, hx, dt, nstep, x);
printMagnetization(l,x,std::cout);
}
return 0;
}
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment