Commit 9443640b authored by Florian Goth's avatar Florian Goth
Browse files

Merge branch '76-ssh-model' of git.physik.uni-wuerzburg.de:marqov/marqov into 76-ssh-model

parents 5fc4b332 5dd0276a
Pipeline #10278 failed with stage
in 67 minutes and 58 seconds
......@@ -7,48 +7,17 @@ Hamiltonian = SSH
L = 4
Ltime = 20
dim = 2
dim = 1
rep = 1
betaMC = 1.0
betaQM = 1.0
dtau = 0.2
k = 2
m = 1
k = 10
g = 0.0, 0.5, 1, 2
mu = 0.4
[BlumeCapelBipartite]
beta = 0.4, 0.45, 0.50, 0.55, 0.60, 0.65, 0.70, 0.75, 0.80
J = -1.0
DA = 1.0
DB = -1.0
dim = 2
rep = 50
L = 8,12,16,24,32,48,64,96,128
[Phi4]
L = 4,6,8,12,16
rep = 5
dim = 3
beta = 0.678, 0.681, 0.684, 0.687, 0.690, 0.693
lambda = 4.5
mass = 1.0
[Ising]
L = 16
rep = 1
dim = 2
beta = 0.430
J = -1.0
mu = 0.0
g = 0.0
......
......@@ -56,9 +56,8 @@
std::string str_g = "g"+std::to_string(std::get<3>(hp));
std::string str_mu = "mu"+std::to_string(std::get<4>(hp));
double betaQM = std::get<5>(hp);
double dtau = std::get<5>(hp);
int Ltime = std::get<6>(hp);
double dtau = betaQM / double(Ltime);
std::string str_Ltime = "Ltime"+std::to_string(Ltime);
std::string str_dtau = "dtau"+std::to_string(dtau);
......
......@@ -72,6 +72,7 @@ class SSHLattice
retval.push_back(j);
for(int j = i+1; j < this->nsites; ++j)
retval.push_back(j);
break;
}
......
......@@ -9,6 +9,9 @@
#include "../hamparts.h"
#include "../metropolis.h"
#define CREATE_CHI_TABLE // writes out the susceptibility and exits the simulation
//#define USE_CHI_TABLE // read susceptibility from file and use it
// ----------------------------------- OBSERVABLES --------------------------------
......@@ -68,7 +71,7 @@ class SSHTwoPointCorrSpace
const int Ls = grid.len;
const int Lt = grid.lentime;
const double norml = 1.0 / Ls / Ls;
const double norml = 1.0 / Ls / Ls / 2; // division by two accounts for double counting
// std::vector<double> retval;
// for (int i=0; i<10; i++) retval.push_back(0);
......@@ -150,56 +153,76 @@ class SSH_multisite
{
public:
double k, beta, dtau, g, mu;
int L;
int ntau;
int L, ntau;
double *const gdat;
std::complex<double>* dexpk;
double *const ftexp;
SSH_multisite(double g, double mu, double b, double d, int myL) : k(-0.5*g*g), beta(b), dtau(d), L(myL), ntau(std::round(beta/dtau)), g(g), mu(mu),
#ifndef SSH_2D
gdat(new double[ntau*L]),
ftexp(new double[L*L])
#else
gdat(new double[ntau*L*L]),
ftexp(new double[L*L*L*L])
#endif
double *const ftexp;
std::complex<double>* dexpk;
std::vector<double> chi_table;
SSH_multisite(double g, double mu, double b, double d, int myL) : k(-0.5*g*g*d*d),
beta(b),
dtau(d),
L(myL),
ntau(std::round(b/d)),
g(g),
mu(mu),
#ifndef SSH_2D
gdat(new double[ntau*L]),
ftexp(new double[L*L])
#else
gdat(new double[ntau*L*L]),
ftexp(new double[L*L*L*L])
#endif
{
// ntau should be nothing else than Ltime
dexpk = new std::complex<double>[L];
for(int d = 0; d < L; ++d)
{
dexpk[d] = std::exp(std::complex<double>(0.0, -2.0*M_PI/L*d ));
for(int k = 0; k < L; ++k)
ftexp[d*L + k] = std::real(std::exp(std::complex<double>(0.0, -2.0*M_PI*k/L*d )));
}
#ifndef SSH_2D
#ifdef USE_CHI_TABLE
// read susceptibility from file
double innumber;
std::ifstream infile("../log/chi.dat", std::fstream::in);
while (infile >> innumber)
{
chi_table.push_back(innumber);
}
infile.close();
#endif
dexpk = new std::complex<double>[L];
for(int d = 0; d < L; ++d)
{
dexpk[d] = std::exp(std::complex<double>(0.0, -2.0*M_PI/L*d ));
for(int k = 0; k < L; ++k)
ftexp[d*L + k] = std::real(std::exp(std::complex<double>(0.0, -2.0*M_PI*k/L*d )));
}
#ifndef SSH_2D
for(int j = 0; j < L; ++j)
{
double k = (j*2)*M_PI/double(L);
// 1D disperson relation
double eps = -2.0*std::cos(k)-mu;
double k = (j*2)*M_PI/double(L); // Fourier faktor
double eps = -2.0*std::cos(k)-mu; // 1D disperson relation (chain)
for(int dt = 0; dt < ntau; ++dt)
{
gdat[dt * L + j] = 0.5*std::exp((-beta/2 + dt*dtau)*eps)/std::cosh(beta*eps/2);
}
}
#else
#else
for(int jx = 0; jx < L; ++jx)
{
double kx = (jx*2)*M_PI/double(L);
for(int jy = 0; jy < L; ++jy)
{
double ky = (jy*2)*M_PI/double(L);
// 2D disperson relation
double eps = -2.0*std::cos(kx) - std::cos(ky) - mu;
for(int dt = 0; dt < ntau; ++dt)
{
gdat[dt * L*L + L*jx + jy] = 0.5*std::exp((-beta/2 + dt*dtau)*eps)/std::cosh(beta*eps/2);
}
}
for(int jy = 0; jy < L; ++jy)
{
double ky = (jy*2)*M_PI/double(L); // Fourier faktor
double eps = -2.0*std::cos(kx) - 2.0*std::cos(ky) - mu; // 2D disperson relation (square lattice)
for(int dt = 0; dt < ntau; ++dt)
{
gdat[dt * L*L + L*jx + jy] = 0.5*std::exp((-beta/2 + dt*dtau)*eps)/std::cosh(beta*eps/2);
}
}
}
#endif
#endif
}
......@@ -215,22 +238,121 @@ class SSH_multisite
Lattice& grid)
{
double retval = 0;
for (int i=0; i<nbrs.size(); i++)
// write susceptibiliy to file
#ifdef CREATE_CHI_TABLE
cout << "writing susceptibiliy to file ... ";
std::ofstream os("../log/chi.dat");
os << std::fixed << std::setprecision(12);
cout << grid.size() << endl;
for (int i=0; i<grid.size(); i++)
{
auto b1 = svnew-svold;
auto b2 = suscept(grid, rsite, nbrs[i]);
auto b3 = s[i];
for (int j=0; j<grid.size(); j++)
{
os << suscept(grid,i,j) << endl;
}
}
auto a1 = dot(b1, mult(b2,b3));
auto a2 = dot(b3, mult(b2,b1));
// transform to Flo's Mathematica index ordering
//
/*
retval = retval + a1 + a2;
os << std::fixed << std::setprecision(12);
int counter = 0;
int len = grid.len;
int vspace = 2*len*len;
int lentime = grid.lentime;
int dim = 2;
for (int xy1=0; xy1<=len; xy1+=len)
{
for (int i1=xy1; i1<vspace; i1+=len*dim)
{
for (int j1=0; j1<len; j1++)
{
for (int k1=0; k1<lentime; k1++)
{
int idx1 = i1 + j1 + 2*k1*len*len;
counter++;
for (int xy2=0; xy2<=len; xy2+=len)
{
for (int i2=xy2; i2<vspace; i2+=len*dim)
{
for (int j2=0; j2<len; j2++)
{
for (int k2=0; k2<lentime; k2++)
{
int idx2 = i2 + j2 + 2*k2*len*len;
os << suscept(grid,idx1,idx2) << endl;
}
}
}
}
}
}
}
}
*/
os.close();
cout << "done! exiting ... " << endl;
exit(0);
/* used for debugging
std::ofstream oss("../log/green.dat");
retval += dot(svnew, mult(suscept(grid, rsite, rsite), svnew));
retval -= dot(svold, mult(suscept(grid, rsite, rsite), svold));
for (int idx1=0; idx1<2*L*L; idx1++)
{
for (int idx2=0; idx2<2*L*L; idx2++)
{
auto w1 = grid.getcrds(idx1);
auto w2 = grid.getcrds(idx2);
auto c1 = wick<decltype(w1)>({w1[0],w1[2],0}, {w1[1],w1[3],0}, {w2[0],w2[2],0}, {w2[1],w2[3],0});
if (c1 < 1e-14) c1 = 0;
oss << c1 << "\t";
}
oss << endl;
}
oss.close();
exit(0);
*/
#endif
double retval = 0;
for (int i=0; i<nbrs.size(); i++)
{
auto chi = suscept(grid, rsite, nbrs[i]);
auto new1 = dot(svnew, mult(chi,s[nbrs[i]]));
auto old1 = dot(svold, mult(chi,s[nbrs[i]]));
retval = retval + 2*(new1 - old1); // since matrix is symmetric
// retval = retval + (new1 - old1); // since matrix is symmetric
}
/*
for (int i=0; i<nbrs.size(); i++)
{
auto chi = suscept(grid, nbrs[i], rsite);
auto new2 = dot(s[nbrs[i]], mult(chi,svnew));
auto old2 = dot(s[nbrs[i]], mult(chi,svold));
retval = retval + (new2 - old2);
}
*/
auto chi = suscept(grid, rsite, rsite);
auto newself = dot(svnew, mult(chi, svnew));
auto oldself = dot(svold, mult(chi, svold));
retval = retval + (newself-oldself);
return retval;
}
......@@ -261,7 +383,7 @@ class SSH_multisite
// Green's function for 2+1 dimensions
// takes two coordinate vectors (x,y,t)
template <typename VertexType>
[[gnu::hot, gnu::optimize("fast-math") ]]
[[gnu::hot, gnu::optimize("fast-math") ]]
double green(const VertexType& c1, const VertexType& c2)
{
double retval = 0;
......@@ -297,14 +419,19 @@ class SSH_multisite
{
// dispersion relation
// do the summation
int index = dti*L*L + jx*L + jy;
double fourierterm = expk.real() * gdat[dti*L*L + jx*L + jy];
retval += expk.real() * gdat[dti*L*L + jx*L + jy];
// increment Fourier transform
expk *= dexpky;
}
expk *= dexpkx; // increment Fourier transform
}
const double norml = 1.0 / pow(2*L,2); // the number of sites per time slice
const double norml = 1.0 / pow(L,2); // the number of sites (not bonds!)
return norml*signum*retval;
}
......@@ -314,9 +441,10 @@ class SSH_multisite
// Green's function for 1+1 dimensions
// takes two coordinate vectors (x,t)
template <typename VertexType>
[[gnu::hot, gnu::optimize("fast-math"), gnu::pure ]]
[[gnu::hot, gnu::optimize("fast-math"), gnu::pure ]]
inline double green(const VertexType& c1, const VertexType& c2)
{
// space
auto dist = std::lrint(c1[0]-c2[0]);
if (dist < 0) dist = L + dist;
......@@ -341,12 +469,12 @@ class SSH_multisite
// compute Greens function in Fourier space
const double *const __restrict__ dat = ftexp + dist*L;
const double *const __restrict__ gdatloc = gdat + dti*L;
double retval = 0;
const double *const __restrict__ dat = ftexp + dist*L;
const double *const __restrict__ gdatloc = gdat + dti*L;
double retval = 0;
for (int j = 0; j < L; ++j)
{
retval += dat[j] * gdat[j];
retval += dat[j] * gdatloc[j];
}
const double norml = 1.0/L;
......@@ -377,6 +505,14 @@ class SSH_multisite
// the actual susceptibility
// takes two indices, representing two bonds in the system
#ifdef USE_CHI_TABLE
template <class Lattice>
double suscept(Lattice& grid, int idx1, int idx2)
{
return chi_table[idx1*grid.size()+idx2];
}
#else
template <class Lattice>
double suscept(Lattice& grid, int idx1, int idx2)
{
......@@ -389,28 +525,25 @@ class SSH_multisite
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<decltype(w1)>({w1[0],w1[2],t1}, {w1[1],w1[3],t1}, {w2[0],w2[2],t2}, {w2[1],w2[3],t2});
const auto c2 = wick<decltype(w1)>({w1[0],w1[2],t1}, {w1[1],w1[3],t1}, {w2[1],w2[3],t2}, {w2[0],w2[2],t2});
const auto c3 = wick<decltype(w1)>({w1[1],w1[3],t1}, {w1[0],w1[2],t1}, {w2[0],w2[2],t2}, {w2[1],w2[3],t2});
const auto c4 = wick<decltype(w1)>({w1[1],w1[3],t1}, {w1[0],w1[2],t1}, {w2[1],w2[3],t2}, {w2[0],w2[2],t2});
// 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]});
const std::complex<double> K1 = green<decltype(w1)>({w1[0],w1[2],t1},{w1[1],w1[3],t1})
+ green<decltype(w1)>({w1[1],w1[3],t1},{w1[0],w1[2],t1});
const std::complex<double> K2 = green<decltype(w1)>({w2[0],w2[2],t2},{w2[1],w2[3],t2})
+ green<decltype(w1)>({w2[1],w2[3],t2},{w2[0],w2[2],t2});
const std::complex<double> K1 = green<decltype(w1)>({w1[0],w1[2],t1},{w1[1],w1[3],t1}) + green<decltype(w1)>({w1[1],w1[3],t1},{w1[0],w1[2],t1});
const std::complex<double> K2 = green<decltype(w1)>({w2[0],w2[2],t2},{w2[1],w2[3],t2}) + green<decltype(w1)>({w2[1],w2[3],t2},{w2[0],w2[2],t2});
#else
const auto t1 = w1[2];
const auto t2 = w2[2];
std::array<double, 2> v10 = {w1[0], t1};
std::array<double, 2> v11 = {w1[1], t1};
std::array<double, 2> v20 = {w2[0], t2};
std::array<double, 2> v21 = {w2[1], t2};
std::array<double, 2> v10 = {w1[0], t1};
std::array<double, 2> v11 = {w1[1], t1};
std::array<double, 2> v20 = {w2[0], t2};
std::array<double, 2> v21 = {w2[1], t2};
const auto c1 = wick(v10, v11, v20, v21);
const auto c2 = wick(v10, v11, v21, v20);
......@@ -423,6 +556,7 @@ class SSH_multisite
return std::real(c1+c2+c3+c4-K1*K2);
}
#endif
};
......@@ -446,7 +580,7 @@ class SSH_Initializer
StateVector newsv(const StateVector& svold)
{
StateVector retval(svold);
double amp = 0.7;
double amp = 0.5;
double diff = rng.real(-amp, amp);
retval[0] += diff;
return retval;
......@@ -457,7 +591,7 @@ class SSH_Initializer
// ------------------------------ HAMILTONIAN ---------------------------
template <class Lattice, typename SpinType = double>
template <typename SpinType = double>
class SSH
{
public:
......@@ -468,13 +602,18 @@ class SSH
template <typename RNG>
using MetroInitializer = SSH_Initializer<StateVector, RNG>;
Lattice& grid;
static constexpr uint Nalpha = 1;
static constexpr uint Nbeta = 1;
static constexpr uint Ngamma = 1;
SSH(double m, double k, double g, double mu, double bQ, int Ltime, int L, Lattice& lattice) : m(m), k(k), g(g), mu(mu), dtau(bQ/double(Ltime)), L(L), betaQM(bQ), name("SSH"), grid(lattice)
SSH(double m, double k, double g, double mu, double dtau, int Ltime, int L) : m(m),
k(k),
g(g),
mu(mu),
dtau(dtau),
betaQM(dtau*Ltime),
L(L),
name("SSH")
{
interactions[0] = new SSH_interaction<StateVector>(m, dtau);
onsite[0] = new SSH_onsite<StateVector>(m, k, dtau);
......@@ -489,8 +628,9 @@ class SSH
// instantiate and choose observables
SSHMag obs_m;
SSHMagSq obs_msq;
SSHTwoPointCorrSpace obs_corr;
auto getobs() {return std::make_tuple(obs_m, obs_msq, obs_corr);}
// SSHTwoPointCorrSpace obs_corr;
auto getobs() {return std::make_tuple(obs_m, obs_msq);}
// auto getobs() {return std::make_tuple(obs_m, obs_msq, obs_corr);}
// initialize state space
......@@ -499,7 +639,7 @@ class SSH
{
for (int i=0; i<grid.size(); i++)
{
if (rng.real() > 0.1) statespace[i][0] = 1;
if (rng.real() > 0.0) statespace[i][0] = 1;
else statespace[i][0] = -1;
}
}
......
......@@ -82,7 +82,7 @@ void selectsim(RegistryDB& registry, std::string outbasedir, std::string logbase
// ----------------- select simulation ------------------
auto beta = registry.Get<std::vector<double> >("mc.ini", ham, "betaMC");
auto betaQM = registry.Get<std::vector<double> >("mc.ini", ham, "betaQM");
auto dtau = registry.Get<std::vector<double> >("mc.ini", ham, "dtau");
auto m = registry.Get<std::vector<double> >("mc.ini", ham, "m");
auto k = registry.Get<std::vector<double> >("mc.ini", ham, "k");
auto g = registry.Get<std::vector<double> >("mc.ini", ham, "g");
......@@ -98,7 +98,7 @@ void selectsim(RegistryDB& registry, std::string outbasedir, std::string logbase
// we need "L" and "Ltime" as explicit parameters in the Hamiltonian
// which requires some gymnastics ...
std::vector<int> dummy = {0};
auto hp = cart_prod(beta, m, k, g, mu, betaQM, dummy, dummy);
auto hp = cart_prod(beta, m, k, g, mu, dtau, dummy, dummy);
// prepare lattices
......@@ -117,7 +117,7 @@ void selectsim(RegistryDB& registry, std::string outbasedir, std::string logbase
}
latts.emplace_back(lat);
}
typedef decltype(sshfilter(latts[0][0], std::declval<std::pair<MARQOV::Config, decltype(hp[0])>>()
)) PPType;
typedef decltype( finalize_parameter_pair(std::declval<MARQOV::Config>(),
......@@ -126,6 +126,7 @@ void selectsim(RegistryDB& registry, std::string outbasedir, std::string logbase
typename GetSchedulerType<
SSH<SSHLattice, double>, SSHLattice, PPType>::MarqovScheduler sched(1);
for (std::size_t j=0; j<nL.size(); j++)
{
for (std::size_t jj=0; jj<nLtime.size(); jj++)
......@@ -142,16 +143,17 @@ void selectsim(RegistryDB& registry, std::string outbasedir, std::string logbase
std::string outpath = outbasedir+"/"+std::to_string(L)+"/";
MARQOV::Config mp(outpath);
mp.setnsweeps(1);
mp.setnsweeps(10);
mp.setncluster(0);
mp.setwarmupsteps(0);
mp.setgameloopsteps(10);
mp.setgameloopsteps(300);
makeDir(mp.outpath);
// set up parameters
auto params = finalize_parameter_pair(mp, hp);
auto rparams = finalize_parameter_pair(mp, hp);
auto params = replicator_pair(rparams, nreplicas[0]);
SSHLattice& latt = latts[j][jj];
auto f = [&latt, &outbasedir, L](auto p){return sshfilter(latt, p);}; //partially apply filter
......@@ -183,11 +185,11 @@ int main()
std::string command;
command = "rm -r " + outbasedir;
system(command.c_str());
command = "rm -r " + logbasedir;
system(command.c_str());
// command = "rm -r " + logbasedir;
// system(command.c_str());
makeDir(outbasedir);
makeDir(logbasedir);
// makeDir(logbasedir);
selectsim(registry, outbasedir, logbasedir);
}
......@@ -92,6 +92,7 @@ namespace MARQOV
template <class StateSpace, class M, class RNGType>
int Metropolis<Hamiltonian, Lattice>::move(const Hamiltonian& ham, const Lattice& grid, StateSpace& statespace, M& metro, RNGCache<RNGType>& rng, double beta, int rsite)
{
typedef typename Hamiltonian::StateVector StateVector;
// old state vector at rsite
......@@ -158,7 +159,6 @@ namespace MARQOV
double dE = interactionenergydiff + onsiteenergydiff + multisiteenergdiff;
int retval = 0;
if ( dE <= 0 )
{
......