SM_MonteCarloSystem.hpp

#ifndef SM_MonteCarloSystem_HPP
#define SM_MonteCarloSystem_HPP
 
template<class StochOutputImplement>
inline tBoolean SM_MonteCarloSystem::stochasticRun(const tIndex& steppersNumber,
                                                   SM_MultiStochasticFunctionsInterface& sfs,
                                                   SM_StochasticOutput<StochOutputImplement>& statistics) {
 
      
    //number of drawn for stochastic output 
    const tInteger &nDrawnSteps=statistics.getDrawnStepsNumber();
 
    //no statictics
    if (nDrawnSteps==0) return stochasticRun(steppersNumber,sfs);
 
    tString routineName="SMOMPI_MonteCarloSystem::stochasticRun<StochOutputI>(const tIndex&,SM_StochasticFunction&,StochOutputI&)";
    CORE_Profiler::StartCallRoutine(routineName);
 
        
    //return value
    tBoolean succeeds=true;
        
    //network of the system
    const SM_Network& network=this->getNetwork();
        
    //material of the system
    const SM_Material& material=this->getMaterial();
        
    //number of particles
    const tIndex& nParticles=network.getParticlesNumber();
    //number of halo particles
    const tIndex& nHaloParticles=network.getHaloParticlesNumber();
        
    //iterator atomic spin direrctions
    SM_RealField& St=getMagneticMomentDirections();
    if (St.getElementsNumber()!=nParticles+nHaloParticles) {
        throw CORE_Exception("stochmagnet/montecarlo",
                             "SMOMPI_MonteCarloSystem::stochasticRun<StochF>()",
                             "bad size of magnetic moments directions");
    }
    //energy scaling  alpha= \frac{\bar E }{k_B.T} 
    const tReal& alpha=getEnergyScaling();
        
    //amplitude of the trial move
    //sigma=\frac{2}{25}*\left( \frac{T*k_B}{\mu_B} \right))^{\frac{1}{5}}
    const tReal& sigma=getNoiseRate();
        
    //get the step index
    tIndex& t=getStepIndex();
    tIndex iStep=0;
    if (nDrawnSteps==1) {
        while (iStep<steppersNumber) {
                
            MCStep(sfs,material,network,St,t,
                   alpha,sigma,mMovesNumber,mRejectsNumber);
                
            //store the statictics
            statistics.store(*this);
            //next step  
            updateStateForNextStep(St);
            t++;
            iStep++;
        }//end loop on stepper of equilibrium      
    }
    else {
        tInteger i;
            
        do {
            i=0;
            while (i<nDrawnSteps) {
                MCStep(sfs,material,network,St,t,
                       alpha,sigma,mMovesNumber,mRejectsNumber);
                //next step
                updateStateForNextStep(St);
                i++;
                t++;
                iStep++;
            }
            //store the statistics
            statistics.store(*this);
        } while (iStep<steppersNumber);
    }
      
    CORE_Profiler::EndCallRoutine(routineName);
    return succeeds;
}
    
template<class StochFunctionsImplement>
inline tBoolean SM_MonteCarloSystem::stochasticRunWithSCStochasticFunctions(const tIndex& steppersNumber,
                                                                            SM_MultiStochasticFunctions<StochFunctionsImplement>& sfs) {
    tString routineName="SMOMPI_MonteCarloSystem::templatedStochasticRun<StochFI>(const tIndex&,Multi<StochF>&)";
    CORE_Profiler::StartCallRoutine(routineName);
        
    //return value
    tBoolean succeeds=true;
        
    //network of the system
    const SM_Network& network=this->getNetwork();
        
    //material of the system
    const SM_Material& material=this->getMaterial();
        
    //number of particles
    const tIndex& nParticles=network.getParticlesNumber();
 
    //number of halo particles
    const tIndex& nHaloParticles=network.getHaloParticlesNumber();
        
    //iterator atomic spin direrctions
    SM_RealField& St=getMagneticMomentDirections();
    if (St.getElementsNumber()!=nParticles+nHaloParticles) {
        throw CORE_Exception("stochmagnet/montecarlo",
                             "SMOMPI_MonteCarloSystem::templatedStochasticRun<StochF>()",
                             "bad size of magnetic moments directions");
    }
      
        
    //energy scaling  alpha= \frac{\bar E }{k_B.T} 
    const tReal& alpha=getEnergyScaling();
        
    //amplitude of the trial move
    //sigma=\frac{2}{25}*\left( \frac{T*k_B}{\mu_B} \right))^{\frac{1}{5}}
    const tReal& sigma=getNoiseRate();
        
    //get the step index
    tIndex& t=getStepIndex();
    tIndex iStep=0;
    while (iStep<steppersNumber) {
            
        MCStepWithSCStochasticFunctions<StochFunctionsImplement>(sfs,material,network,St,t,
                                                                 alpha,sigma,mMovesNumber,mRejectsNumber);
        //next step
        updateStateForNextStep(St);
        t++;
        iStep++;
    }//end loop on stepper of equilibrium      
        
    CORE_Profiler::EndCallRoutine(routineName);
    return succeeds;
        
}
    
template<class StochFunctionsImplement,class StochOutputImplement>
inline tBoolean SM_MonteCarloSystem::stochasticRunWithSCStochasticFunctions(const tIndex& steppersNumber,
                                                                            SM_MultiStochasticFunctions<StochFunctionsImplement>& sfs,
                                                                            SM_StochasticOutput<StochOutputImplement>& statistics) {
        
        
        
    //number of drawn for stochastic output 
    const tInteger &nDrawnSteps=statistics.getDrawnStepsNumber();
 
    //no statictics
    if (nDrawnSteps==0) return stochasticRunWithSCStochasticFunctions<StochFunctionsImplement>(steppersNumber,sfs);
 
    tString routineName="SMOMPI_MonteCarloSystem::templatedStochasticRun<StochFI,StochOI>(const tIndex&,StochF&,StochOI&)";
    CORE_Profiler::StartCallRoutine(routineName);
        
    //return value
    tBoolean succeeds=true;
        
    //network of the system
    const SM_Network& network=this->getNetwork();
        
    //material of the system
    const SM_Material& material=this->getMaterial();
        
    //number of particles
    const tIndex& nParticles=network.getParticlesNumber();
 
    //number of halo particles
    const tIndex& nHaloParticles=network.getHaloParticlesNumber();
        
    //iterator atomic spin direrctions
    SM_RealField& St=getMagneticMomentDirections();
    if (St.getElementsNumber()!=nParticles+nHaloParticles) {
        throw CORE_Exception("stochmagnet/montecarlo",
                             "SMOMPI_MonteCarloSystem::templatedStochasticRun<StochF,StochO>()",
                             "bad size of magnetic moments directions");
    }
    //energy scaling  alpha= \frac{\bar E }{k_B.T} 
    const tReal& alpha=getEnergyScaling();
        
    //amplitude of the trial move
    //sigma=\frac{2}{25}*\left( \frac{T*k_B}{\mu_B} \right))^{\frac{1}{5}}
    const tReal& sigma=getNoiseRate();
        
    //get the step index
    tIndex& t=getStepIndex();
    tIndex iStep=0;
    if (nDrawnSteps==1) {
        while (iStep<steppersNumber) {
                
            MCStepWithSCStochasticFunctions<StochFunctionsImplement>(sfs,material,network,St,t,
                                                                     alpha,sigma,mMovesNumber,mRejectsNumber);
                
            //store the statictics
            statistics.store(*this);
            //next step
            updateStateForNextStep(St);
            t++;
            iStep++;
        }//end loop on stepper of equilibrium      
    }
    else {
        tInteger i;
        do {
            i=0;
            while (i<nDrawnSteps) {
                MCStepWithSCStochasticFunctions<StochFunctionsImplement>(sfs,material,network,St,t,
                                                                         alpha,sigma,mMovesNumber,mRejectsNumber);
 
                //next step
                updateStateForNextStep(St);
                i++;
                t++;
                iStep++;
            }
            //store the statistics
            statistics.store(*this);
        } while (iStep<steppersNumber);
    }
    CORE_Profiler::EndCallRoutine(routineName);
    return succeeds;
}
 
 
template<class StochFunctionsImplement>
inline void SM_MonteCarloSystem::MCStepWithSCStochasticFunctions(SM_MultiStochasticFunctions<StochFunctionsImplement>& multiStochasticFunctions,
                                                                 const SM_Material& material,const SM_Network& network,
                                                                 SM_RealField& S,
                                                                 const tIndex& stepIndex,const tReal& alpha,const tReal& sigma,
                                                                 tIndex& nMoves,tIndex& nRejects)  {
        
   
    //number of particles, ignore halo particles
    const tIndex& nParticles=network.getParticlesNumber();
    
    tReal u=0;//uniform random in [0,1[
    tIndex i=0;//index of the spin
    tReal *Si=null;//atomic spin direction of the spin i
    const tReal *eSi=null;//end ietartor at coordinate of spin
    std::array<tReal,SM_Constants::DIM> oldSi;//od values of Si  before moving
    tReal Eold;//old energy of the spin i before moving
       
    tReal P;//probability to accepet the move
    //stochastic function for threadId 0
    SM_StochasticFunctions<StochFunctionsImplement>& stochF=multiStochasticFunctions[0];
        
    for(tIndex r=0;r<nParticles;r++) {
            
        //number of moves
        nMoves++;
            
        //randomize the spin
        u=stochF.scUniformRandom();
            
        //index of the spin
        i=nParticles*u;
            
        //direction of the spin i
        Si=&S[0];
        Si+=SM_Constants::DIM*i;
        eSi=Si;eSi+=SM_Constants::DIM;
            
        //store the S value at i into oldS
        memcpy(oldSi.data(),Si,SM_Constants::DIM*sizeof(tReal));  //dest,source,size
 
        //update state of operators
        this->updateOperatorsState(stepIndex,network,material,S);
            
        //compute old spin energy
        Eold=this->computeSpinEnergy(i,stepIndex,network,material,S);
            
        //trial move of the spin
        this->randomTrialMoveWithSCStochasticFunctions(stochF,sigma,Si,eSi);
            
        //compute energy of new state of the spin i
        P=this->computeSpinEnergy(i,stepIndex,network,material,S);
            
        //compute probability to keep the new change
        //P=exp(-\alpha.\Delta E)
        P-=Eold;
        P*=alpha;
        P*=-1.L;
        P=exp(P);
            
        u=stochF.scUniformRandom();
        //accept only if P>=u
        if (u>P) {
            //reject the new position
            memcpy(Si,oldSi.data(),SM_Constants::DIM*sizeof(tReal));  //dest,source,size
            nRejects++;
        }
            
    }//end loop on particles drawn
}
 
template<class StochFunctionsImplement>
inline void SM_MonteCarloSystem::randomTrialMoveWithSCStochasticFunctions(SM_StochasticFunctions<StochFunctionsImplement>& stochasticFunctions,
                                                                          const tReal& sigma,
                                                                          tReal *bS,const tReal* eS) const{
    tReal *S=bS;
    tReal invNormS=0;
        
    tFlag moveType=(tFlag)(3.*stochasticFunctions.scUniformRandom());//type of move
        
    switch(moveType) {
            
    case 0://flip move
        while (S!=eS) {
            (*S)*=-1;
            S++;
        }
        break;
    case 1://gaussian move
            
        while (S!=eS) {
            (*S)+=sigma*stochasticFunctions.scNormalRandom();
            invNormS+=(*S)*(*S);
            S++;
        }
        invNormS=1./sqrt(invNormS);
        S=bS;
        while (S!=eS) {
            (*S)*=invNormS;
            S++;
        }
        break;
    case 2://random move
            
        while (S!=eS) {
            (*S)=stochasticFunctions.scNormalRandom();
            invNormS+=(*S)*(*S);
            S++;
        }
        invNormS=1./sqrt(invNormS);
        S=bS;
        while (S!=eS) {
            (*S)*=invNormS;
            S++;
        }
        break;
            
    }
}
  
#endif