Hi Daniel,

Indeed nobody ever wanted to change g and some modification dropped the energy dependence on g. The following patch should fix it.

Matthias

Index: localupdate.h =================================================================== --- localupdate.h (revision 6036) +++ localupdate.h (working copy) @@ -42,7 +42,7 @@ template <class Graph, class MomentMap, class RNG, class CouplingMap, class SpinfactorMap> void local_update(const Graph& graph, MomentMap& moment, double beta, RNG& rng, const CouplingMap& coupling, - const SpinfactorMap& spinfactor, + const SpinfactorMap& spinfactor, double g, TinyVector<double,CouplingMap::value_type::dim> h = TinyVector<double,CouplingMap::value_type::dim>(0.0))

@@ -72,7 +72,7 @@ delta_E += bond_energy_change(moment[boost::source(*edge,graph)], moment[boost::target(*edge,graph)],J,update); } - delta_E += spinfactor[vertex]*site_energy_change(moment[vertex], h, update); + delta_E += spinfactor[vertex]*site_energy_change(moment[vertex], h, update) * g;

// step iv) if (delta_E<0 || rng() < exp(-beta*delta_E)) @@ -86,7 +86,7 @@ void local_update_self(const Graph& graph, MomentMap& moment, double beta, RNG& rng, const CouplingMap& coupling, const SelfinteractionMap& selfinteraction, - const SpinfactorMap& spinfactor, + const SpinfactorMap& spinfactor, double g, TinyVector<double,CouplingMap::value_type::dim> h = TinyVector<double,CouplingMap::value_type::dim>(0.0)) { @@ -115,7 +115,7 @@ delta_E += bond_energy_change(moment[boost::source(*edge,graph)], moment[boost::target(*edge,graph)],J,update); } - delta_E += spinfactor[vertex]*site_energy_change(moment[vertex], h, update); + delta_E += spinfactor[vertex]*site_energy_change(moment[vertex], h, update) * g;

delta_E += onsite_energy_change(moment[vertex], selfinteraction[vertex],update);

Index: abstractspinsim.h =================================================================== --- abstractspinsim.h (revision 6036) +++ abstractspinsim.h (working copy) @@ -205,7 +205,7 @@ MAT this_J; this_J.setMatrix(J_name,parms); bond_type[*b] = num_bond_types; - this_J = this_J /(g_*g_) * (quantum_convention_ ? 1. : -1.); + this_J = this_J * (quantum_convention_ ? 1. : -1.); couplings_[*b]=this_J; double det = this_J.det(); couplings_det[*b] = det; Index: spinsim.h =================================================================== --- spinsim.h (revision 6036) +++ spinsim.h (working copy) @@ -166,10 +166,10 @@ typedef AbstractSpinSim<MAT> parent; if (this->general_case_) local_update_self(this->graph(),spins_, this->beta_,this->random_01, - this->couplings_,this->selfinteraction_,this->spinfactor_,this->h_); + this->couplings_,this->selfinteraction_,this->spinfactor_,this->g_,this->h_); else local_update(this->graph(),spins_, this->beta_,this->random_01, - this->couplings_,this->spinfactor_,this->h_); + this->couplings_,this->spinfactor_,this->g_,this->h_); update_info_type up_info; up_info.clustersize = this->num_sites(); up_info.m2_ = 0; @@ -190,10 +190,10 @@ else { if (this->general_case_) local_update_self(this->graph(),spins_,this->beta_,this->random_01, - this->couplings_,this->selfinteraction_,this->spinfactor_,this->h_); + this->couplings_,this->selfinteraction_,this->spinfactor_,this->g_,this->h_); else local_update(this->graph(),spins_,this->beta_,this->random_01, - this->couplings_,this->spinfactor_,this->h_); + this->couplings_,this->spinfactor_,this->g_,this->h_); update_info_type up_info; up_info.clustersize = this->num_sites(); up_info.m2_ = 0;

On 7 Mar 2012, at 07:43, daniel.aravena@antares.qi.ub.edu wrote:

Hi Matthias,

I tried the new patch and it works correctly, the only thing I find strange is that the "Magnetization along field" vs. magnetic field plot doesn't look like dependent on the g factor but the energies are dependent on g, so i would expect a higher magnetization for higher g at the same field when far from saturation while i get the same magnetization regardless of the value of g (i tried g=1.0 and 2.0 and h from 0 to 11).

Thanks Kind Regards Daniel

...

Hi Daniel,

Can you please try to apply the following patch, and turn the quantum convention off. It will not give sensible results in magnetic fields.

Matthias

Index: applications/mc/spins/abstractspinsim.h

--- applications/mc/spins/abstractspinsim.h (revision 6011) +++ applications/mc/spins/abstractspinsim.h (working copy) @@ -136,6 +136,11 @@ else use_error_limit = false;

• if (has_magnetic_field_ && quantum_convention_)
• boost::throw_exception(

•  std::runtime_error("Cannot use quantum convention with nonzero
  

magnetic field ."));

if (parms["UPDATE"]=="cluster" && has_magnetic_field_) {boost::throw_exception(std::runtime_error("Illegal update type " + std::string(parms["UPDATE"]) @@ -159,7 +164,7 @@ spinvalue = quantum_convention_ ? 0.5 : 1.; else spinvalue = alps::evaluate<double>(parms[spin],parms);

spinfactor_[*si]=std::sqrt(spinvalue*(spinvalue+(quantum_convention_ ? 1. : 0.)))*g_;

spinfactor_[*si]=std::sqrt(spinvalue*(spinvalue+(quantum_convention_ ? 1. : 0.)));

 // set the values for the diagonal matrix terms
 site_type_visited[type]=true;

@@ -176,7 +181,7 @@ selfinteraction_[*si].setMatrix(diag_string,parms); general_case_=true; }

•  selfinteraction_[*si] =
  

selfinteraction_[*si]*spinvalue*(spinvalue+(quantum_convention_ ? 1. : 0.));

•  selfinteraction_[*si] = selfinteraction_[*si]*spinvalue*spinvalue;
  
  

  if (!selfinteraction_[*si].is_symmetric()) boost::throw_exception(std::runtime_error("Invalid matrix for

D")); Index: applications/mc/spins/spinsim.h =================================================================== --- applications/mc/spins/spinsim.h (revision 6011) +++ applications/mc/spins/spinsim.h (working copy) @@ -246,14 +246,14 @@ m_h = this->spinfactor_[*si]*spins_[*si].mag_h(this->h_normalized);

if (this->has_magnetic_field_)

• en+=this->spinfactor_[*si]*site_energy(spins_[*si],this->h_);

en+=this->spinfactor_[*si]*site_energy(spins_[*si],this->h_)*this->g_;

double par=this->parity(*si); ++si;

for (; si !=this->sites().second ; ++si) { if (this->has_magnetic_field_)

•  en+=this->spinfactor_[*si]*site_energy(spins_[*si],this->h_);
  

en+=this->spinfactor_[*si]*site_energy(spins_[*si],this->h_)*this->g_;

if (this->general_case_) en+=onsite_energy(spins_[*si],this->selfinteraction_[*si]); @@ -321,7 +321,6 @@ if (this->print_sweeps_) { int good_sweeps = this->sweeps_done_-this->thermalization_sweeps_; if (good_sweeps > 0 && (good_sweeps % this->print_sweeps_ ==0)) {

•  std::cout << "Dumping configuration: \n";
  

  for (int i=0;i<this->num_sites();++i) { std::cout << i << " "; std::vector<double> coord = this->coordinate(i);

@@ -350,7 +349,7 @@ bt_energy.resize(bt+1); bt_count.resize(bt+1); }

• double be =

bond_energy(spins_[this->source(*bi)],spins_[this->target(*bi)],

• double be =

this->spinfactor_[this->source(*bi)]*this->spinfactor_[this->target(*bi)]*bond_energy(spins_[this->source(*bi)],spins_[this->target(*bi)], this->couplings_[*bi]);

en -= be; @@ -374,13 +373,16 @@ mst *= this->spinfactor_[*si]; double m_h = 0.; if (this->has_magnetic_field_)

•  {
   m_h =
  

this->spinfactor_[*si]*spins_[*si].mag_h(this->h_normalized);

en+=this->spinfactor_[*si]*site_energy(spins_[*si],this->h_)*this->g_;

•  }
  

double par=this->parity(*si); ++si;

for (; si !=this->sites().second ; ++si) { if (this->has_magnetic_field_)

•  en+=site_energy(spins_[*si],this->h_);
  

en+=this->spinfactor_[*si]*site_energy(spins_[*si],this->h_)*this->g_; if (this->general_case_) en+=onsite_energy(spins_[*si],this->selfinteraction_[*si]);

On Mar 1, 2012, at 10:39 PM, daniel.aravena wrote:

...

ok, here it is the parameter file for S=1/2 and 2 sites:

LATTICE_LIBRARY=/home/g1daniel/ALPS_LATT/latt02.xml MODEL="Heisenberg" LATTICE="2_graph" UPDATE="local" CONVENTION="quantum" THERMALIZATION=1000000 SWEEPS=10000000 S0=1/2 J0=0 g=1.0 T=1.0 {h=0} {h=1} {h=10} {h=100} {h=1000} {h=10000} {h=100000} {h=1000000} {h=10000000} {h=100000000} {h=1000000000} {h=10000000000}

the last field (h=10000000000) gives -1e+10 for the energy

Kind regards Daniel

On Thu, 01 Mar 2012 22:31:43 +0100, Matthias Troyer wrote:

...

No, I'm sorry - I checked it and we actually do include S in the field term as we should. TO look into this Iw ould like to ask you to send me your complete parameter file for one of the cases, e.g. for

...

...

> Spin Number_of_sites Energy > 1/2 2 -1e+10

On 1 Mar 2012, at 22:14, daniel.aravena wrote:

...

perfect, but, how can you handle the field dependence of a system composed of sites of different spins assigning just a single field

Matthias