A well-known example of classical Monte-Carlo method application is Ising model with $S=1/2$.
As I understood, people there widely use it for any kind of magnetic materials following the same idea $$ P_{I\rightarrow f}\propto e^{-\Delta E_{fi}/kT} $$ right?
What does it mean practically? We have an arbitrary Hamiltonian of a large system with two-particles interaction $$ H = H_0 + \sum_{ij} V_{ij}, \quad V_{ij}=I_{ij}\hat{A}_i \hat{A}_j $$ then, we apply sort of mean field approximation and for one subsystem (spin) one can write $$ H_i \approx H_0^{(i)} + \sum_{j} I_{ij} \langle A_j\rangle \hat{A}_i $$ where $\langle A_j \rangle$ are just numbers stored somewhere from the previous calculations.
Ok then, according to the method, we pick up a random $i$ and try to update it's state. In order to do that we find all eigenvalues and eigenstates of $H_i$, i.e. $$ H_i: \quad set~of~\{E_k^{(i)},\Psi_k^{(i)}\} $$ At this stage I do not know what is the correct way of dealing with it using local update strategy
- Pick up a new state $\Psi_m$ for $i$-subsystem by random choice from the subset of $\{E_k,\Psi_k\}$ with weights $e^{- E_k/kT}$ and store $\langle A_j \rangle = \langle \Psi_m | A | \Psi_m \rangle$ for further calculations
- Perform a transition from some initial (previous) state $\Psi_h$ to any of available states $\{E_k, \Psi_k\}$ with probability $e^{-(E_m-E_h)/kT}$ (where $E_h$ is recalculated using the present Hamiltonian on the old $\Psi_h$). And then store $\langle A_j \rangle = \langle \Psi_m | A | \Psi_m \rangle$ as well as a new state $\Psi_m$ for the further calculations
The second (2) approach looks more sophisticated for me and kinda feels more natural due to the transitions caused by the thermal fluctuations, not clear how to implement it.
Due to the lack of knowledge I followed (1) a long time ago and tried it for the ferrimagnetically bonded cubic spin-system with (L=2, S=2) and (L=0, S=3/2). It reproduced quite well the phase transition, hysteresis loop and etc. But I still have doubts in this regard, if this procedure (1) is statistically correct.
I hope some brainy Monte-Carlo guys here can clarify it a bit for me.
PS: Some extra details
For the simplicity I tested it on a simple cubic spin system with S=3/2 $$ H = H_{anisotropy} + H_{dipole-dipole} + H_{exchange} $$
it forms a multidomain state, similar to Ref, where they managed to form magnetic-domains using Classical Monte-Carlo.
It looks right for now...