Motility-Induced Phase Separation(MIPS)

Grupo: Bernardo Boatini e Murilo Kessler Azambuja

Neste trabalho iremos seguir de perto o artigo "Critical motility-induced phase separation belongs to the Ising universality class" ^[1] de Benjamin Partridge e Chiu Fan Lee.

Introdução

A matéria ativa é um tipo de sistema fora do equilíbrio termodinâmico, onde cada "partícula" ou "agente" do sistema tem a capacidade dissipar energia na forma de forças mecânicas exercidas sobre o ambiente no qual está imerso. Esses sistemas muitas vezes podem exibir vários fenômenos novos, como movimentos coletivos quando exige-se um alinhamento das partículas (como bio-polímeros que se auto organizam, tipo os microtúbulos que são parte do citoesqueleto celular, ou cardumes de peixe, por exemplo), ou os chamados MIPS (Motility-Induced Phase Separation) que são sistemas que apresentam uma mudança de fase física devido à interações entre partículas que proíbem a ocupação simultânea de um volume do espaço além delas apresentarem persistência em uma certa direção de movimento. Nestes sistemas, os agentes possuem a capacidede de auto-propulsão, de forma que a condição de balanço detalhado é quebrada no nível microscópico ^[1], uma vez que as partículas possuem uma direção preferencial de movimento.Estes tipos de sistema muitas vezes podem ser tratados a partir de uma abordagem hidrodinâmica, na qual são chamados de fluidos ativos, e divergem do comportamento usual de fluidos compostos de matéria usual inativa, os quais são descritos pela equação de Navier-Stokes.

Apesar desses sistemas fundamentalmente quebrarem a condição do balanço detalhado, ainda não é claro se o comportamento estatístico universal de sistemas de fluidos ativos necessariamente divergem daqueles de sistemas em equilíbrio. A investigação de comportamentos universais nestes tipos de sistema, além de ser um interece central na física, também pode nos permitir utilizar conhecimentos já bem estabelecidos (como a transição de fase em sistemas de matéria em equilíbrio termodinâmico) para descrever sistemas novos. Neste trabalho foi reproduzido o artigo ^[1] onde foi investigada o comportamento crítico de um sistema com MIPS e se viu que este tipo de sistema pertence à classe de universalidade de sistemas em equilíbrio como o modelo de Ising para spin. Para tal, foi utilizada três abordagens diferentes: um modelo hidrodinâmico, uma descrição baseada em teoria de campos e a simulação de um modelo de rede hexagonal.

Modelo

O sistema geral que estamos interessados em descrever é constituído de um conjunto de partículas em um meio com friccção, i.e, em um sistema sem conservação de momentum e com interações que não permitem a ocupação simultânea de duas partículas em um mesmo sítio (FIG. 1).

FIGURA 1: Imagem retirada de ^[1]. Em (a) vemos a difusão devido à rotação da direção de movimento preferencial da partícula, em (b) temos o "movimento balístico" (seta verde) e a difusão devido à translação da partícula (setas azuis) e, em (c), vemos a interação entre partículas em sítios vizinhos, impedindo a ocupação simultânea do mesmo sítio por duas partículas.

Uma forma de modelar este sistema é considerá-lo um fluido ativo polar e compressível sem interações de alinhamento entre partículas, de forma que não se espera observar movimentos coletivos. Tal sistema é descrito de forma contínua pela chamada equação de Toner-Tu ^[1].

De um ponto de vista microscópico, podemos usar uma formulação discreta baseada na teoria de campos. Neste caso, consideramos uma rede hexagonal bidimensional, com o numero de ocupação de cada sítio da rede restringido por uma constante $M$ (futuramente diremos que $M=1$ , de forma que apenas uma partícula pode ocupar cada sítio).

Modelo no Continuo

As equações de Tuner-Tu para o modelo em questão são dadas por ^[1]

\partial _{t}\rho +\nabla \cdot \mathbf {g} =0,\qquad (1a)

$\partial _{t}\mathbf {g} =-\zeta \nabla \rho -\kappa \mathbf {g} +\mu \nabla ^{2}\mathbf {g} +\mathbf {f} ,\qquad (1b)$

onde $\rho$ é o campo de densidade de massa e $\mathbf {g}$ é o campo de densidade de momentum. A equação (1a) diz respeito à conservação de massa do sistema como um todo. A equação (1b) se assemelha muito à equação de Navier-Stokes, assim podemos ver que o primeiro termo da equação caracteriza a compressibilidade do fluido, o segundo termo representa uma difusão do campo de momentum (o sinal negativo indica que o campo de momentum diminui com o passar do tempo), o terceiro termo representa a viscosidade do fluido e $\mathbf {f}$ é um ruído Gaussiano com estatísticas espaço-temporais dadas por

$\langle \mathbf {f} (t,\mathbf {r} )\rangle =0,\ \langle \mathbf {f} (t,\mathbf {r} )\mathbf {f} (t',\mathbf {r'} )\rangle =2D\delta (t-t')\delta (\mathbf {r} -\mathbf {r'} ),$

onde $D$ é a intensidade do ruído.

No mesmo artigo ^[1] os autores argumentam que, devido à ausência de interações de alinhamento entre partículas, a emergência de movimento coletivo é impossível. Por este motivo, temos que o campo de densidade de momentum deve ir à zero no limite hidrodinâmico, de forma que $\mathbf {g}$ não é uma variável livre e pode ser escrita em função de $\rho$ . Desta forma, ^[1] obtêm que

$\mathbf {g} =-\nabla {\Big [}(a_{1}\phi +a_{2}\phi ^{2}+a_{3}\phi ^{3})-(\nabla ^{2}(b\phi )){\Big ]}+\mathbf {f} ,$

de modo que, ao substituir em (1a), temos

$\partial _{t}\phi =\nabla ^{2}{\frac {\delta H}{\delta \phi }}+\nabla \cdot \mathbf {f} ,\qquad (2)$

onde $\phi (\mathbf {r} )=\rho (\mathbf {r} )-\rho _{0}$ para alguma constante $\rho _{0}$ e $H={\frac {a_{1}}{2}}\phi ^{2}+{\frac {a_{2}}{3}}\phi ^{3}+{\frac {a_{3}}{4}}\phi ^{4}+{\frac {b}{2}}{\big (}\nabla \phi {\big )}^{2}$ é o Hamiltoniano de Ginzburg-Landau. A partir da equação (2), vemos que é possível descrever o sistema a partir da manipulação de dois parâmetros principais: o campo de densidade de partículas $\rho$ (ou $\phi$ , a menos de uma constante) e a intensidade $D$ do ruído $\mathbf {f}$ .

Modelo no Espaço Discreto

Como já foi dito antes, no modelo discreto é utilizado uma rede bidimensional hexagonal com número de ocupação máximo $M$ para cada sítio da rede. Os agentes são divididos em seis tipos distintos de partículas ativas, cada uma se move em uma direção $\theta _{i}$ (as diagonais do hexágono) e somente vai saltar para o sítio adjacente ao longo desta direção, com uma taxa que depende da ocupação do sítio em questão. Além disso, é possível que uma partícula se transforme em outra, com uma taxa $\sigma$ , correspondendo assim a um ruído rotacional de cada partícula. Utilizando a notação utilizada em ^[1], temos que a rede hexagonal é caracterizado pelo conjunto de vetores $\mathbf {R}$ e a orientação ou tipo da partícula ativa pelo índice $i$ . Desta forma o número de partículas do tipo $i$ ocupando o sítio $\mathbf {R}$ da rede é denotado por $A_{\mathbf {R} }^{\theta _{i}}$ . Com o formalismo desenvolvido em ^[2], temos que a teoria de campos nos fornece uma ação dada por

$S=\int dt\sum _{\mathbf {R} }{\Bigg \{}\sum _{i}{\Bigg [}{\hat {A}}_{\mathbf {R} }^{\theta _{i}}\partial _{t}A_{\mathbf {R} }^{\theta _{i}}$

$-\gamma {\Big (}M-N_{\mathbf {R} +\mathbf {e} _{\theta _{i}}}{\Big )}A_{\mathbf {R} }^{\theta _{i}}{\Big (}e^{{\hat {A}}_{\mathbf {R} +\mathbf {e} _{\theta _{i}}}^{\theta _{i}}-{\hat {A}}_{\mathbf {R} }^{\theta _{i}}}-1{\Big )}{\Bigg ]}\qquad (3)$

$-\sum _{\langle i,j\rangle }\sigma A_{\mathbf {R} }^{\theta _{i}}{\Big (}e^{{\hat {A}}_{\mathbf {R} }^{\theta _{j}}-{\hat {A}}_{\mathbf {R} }^{\theta _{i}}}-1{\Big )}{\Bigg \}},$

onde $N_{\mathbf {R} }=\sum _{j}A_{\mathbf {R} }^{\theta _{j}}$ é o número de partículas total ocupando o sítio $\mathbf {R}$ . A segunda linha da equação (3) diz respeito ao termo de transição entre uma posição $\mathbf {R}$ e uma posição $\mathbf {R} +\mathbf {e} _{\theta _{i}}$ na rede hexagonal. Esta transição ocorre à uma taxa $\gamma {\Big (}M-N_{\mathbf {R} +\mathbf {e} _{\theta _{i}}}{\Big )}$ , de forma que a taxa de transição para o sítio decresce com o aumento do número de partículas ocupando o mesmo. Se o número de partículas ocupando o sítio é $M$ , então a taxa de transição é nula, indicando que não será possível que ocorra uma transição para este sítio. Este termo modela as interações de exclusão de volume discutidas anteriormente. Em nosso caso específico tomaremos $M=1$ . Além disso, na terceira linha da equação (3) corresponde à transição de tipo da partícula, i.e., à transição do movimento preferencial da partícula de uma orientação $\theta _{i}$ para uma $\theta _{j}$ . Esta transição ocorre com uma taxa $\sigma$ e a notação $\langle i,j\rangle$ indica que a partícula somente trocará para uma orientação $\theta _{j}$ que seja adjascente à orientação $\theta _{i}$ inicial.

Implementação

Nessa secção serão colocados alguns detalhes da implementação da rede, das medidas e da dinâmica de Monte Carlo.

Algoritmo

A dinâmica de Monte Carlo foi implementada de forma que a cada MCS, N partículas do sistema são selecionadas aleatoriamente:

- Primeiramente a partícula decide se vai mudar seu sentido( $\theta$ ) sorteando um valor $d$ através de uma distribuição gaussiana centrada em 0, com variância igula a intencidade do ruído $\sigma$ ;

- A nova direção é dada por Falhou ao verificar gramática (erro de sintaxe): {\displaystyle \theta + d.60°} ;

- Depois disso a partícula tem probabilidade $24/30$ de andar para o sitio vizinho naquela direção(se ele estiver vazio);

- Caso a partícula não decida andar naquela direção, ela se desloca na direção de outro sitio vizinho(se ele estiver vazio) sorteado aleatoriamente(emulando o comportamento difusivo).

Redes Hexagonais

Em modelos do tipo LGCA(Lattice Gas Celular Automata) o tipo de rede, bem como as topologias usadas, podem ser determinantes para as medidas e caracterizações do estado macroscópico. A simetria hexagonal se mostra uma boa alternativa para aumentar os graus de liberdade de uma partícula numa rede, quando comparada com a simetria quadrada. Na geometria hexagonal, assim como na quadrada, todos os sítios vizinhos possuem uma mesma distancia entre si só que com 6 graus de liberdade ao invés de quatro.

FIGURA 2: Imagem adaptada de ^[3]. As setas duplas em verde sinalizam os sentidos que seguem trajetórias fechadas passando por só uma borda na topologia "armchair"(a) e "zigzag" (b)

Existem duas topologias possíveis em uma rede hexagonal: "zigzag" e "armchair"^[3], e a diferença das duas está principalmente na implementação das condições de contorno periódicas. No caso da topologia em "zigzag"(FIG.2b), o espaço é representado em um paralelogramo, e existem 4 sentidos em que cruzando uma parede a partícula fecha sua trajetória sem ter que cruzar nova mente outra parede do sistema. Na topologia "armchair"(FIG.2a) por outro lado só existem 2 sentidos em que isso acontece.

Sub-Box Sampling Method

O método das subcaixas consiste em dividir o sistema em uma certa quantidade de subsistemas de interesse, permitindo observar as flutuações de densidades, para um valor fixo de partículas no sistema, apenas nessas regiões.

Nesse trabalho o sistema foi dividido em 4 caixas de tamanho $L\times L$ : 2 centradas no centro de massa horizontal do sistema e duas centradas na região horizontal "mais rarefeita", que consequentemente esta a uma distancia de $3L$ do centro de massa.

Quando os parâmetros permitem que haja separação de fase, duas subcaixas ficam posicionadas exatamente sobre o centro da fase liquida e a mesma coisa para a fase gasosa, permitindo uma medição mais precisa das flutuações de densidade sem considerar as regiões de fronteira.

Porém, essas flutuações se tornam especialmente relevantes no ponto crítico, onde o uso de subcaixas se torna muito importante na extração medidas para classificar a transição de fase do modelo em análise.

Medida 1

A primeira medida utilizada para comparar o MIPS com o Ising foi a compressibilidade isotérmica dos subsistemas $\chi _{L}$ definida da seguinte forma

$\chi _{L}={\frac {\langle N^{2}\rangle _{L}-\langle N\rangle _{L}^{2}}{\langle N\rangle _{L}}}$

onde $\langle N\rangle _{L}$ denota a média do número de partículas contidas em um subsistema de tamanho L.

A compressibilidade isotérmica possui um fator de escala com o tamanho do sistema( $\chi _{L}\alpha L^{\gamma /\nu }$ ), que é conhecido para o caso do Ising 2D ( $\gamma /\nu =7/4$ ).

Medida 2

A segunda medida proposta é a derivada do cumulante de quarta ordem ( $U_{4}$ ) com relação a distância adimensional do ponto crítico ( $\tau ={\frac {\sigma -\sigma _{c}}{\sigma }}$ ) calculada no ponto crítico( $\tau =0$ ). O cumulante é definido da seguinte forma

$U_{4}={\frac {\langle \Delta \rho ^{2}\rangle _{L}^{2}}{\langle \Delta \rho ^{4}\rangle _{L}}}$

$\Delta \rho =\rho -\langle \rho \rangle _{L}$

Em ^[1] $\langle .\rangle _{L}$ é descrito como: "[...] an average over the four sub-boxes of finite size L [...]"

O cumulante de quarta ordem possui a propriedade de ser invariante frente a escala do sistema no ponto crítico(FIG.3), e por isso a sua derivada no ponto critico se torna um fator de escala importante para comparar o estado macroscópico do MIPS e do Ising.

FIGURA 3: Imagem retirada de ^[1]. Mostra o comportamento esperado do cumulante de quarta ordem para o modelo de MIPS nas proximidades da criticidade em sistemas de diferentes tamanhos.

Em $\tau =0$ a derivada do cumulante de quarta ordem, com respeito a $\tau$ , possui um fator de escala $\alpha L^{1/\nu }$ , que também pode ser encontrado no modelo de Ising( $\nu =1$ ).

Resultados

Com base nos Diagrama de Fase do sistema(FIG.4) foi possivel extrair os parâmetros característicos dos principais comportamentos do sistema, e testar a simulação nesses diferentes casos. O tempo de relaxação do sistema $t_{relax}=5\times 10^{5}MCS$ (medido a partir do tem pode correlação) foi o mesmo utilizado no artigo. a densidade do sistema foi fixada em $\rho =0.522$ em todos as simulações. Embora o artigo ^[1] não especifique a topologia de rede utilizada, nesta análise foi adotada a topologia "armchair".

FIGURA 4: Imagem retirada de ^[1]. Mostra o diagrama de fase do sistema.

No estado de equilíbrio, após o $t_{relax}$ , para $\sigma =0.34$ tem-se um comportamento homogêneo(FIG.5), onde o ruído é tão intenso que nenhuma fase é formada no sistema.

FIGURA 5:Imagem de simulação por Bernardo Boatini. Sistema Homogêneo no estado estacionário. L=36

Quando $\sigma =0.25$ é possivel simular um estado estacionário, após $t_{relax}$ , em que ocorre separação de fase. Na figura 6 é evidente a distinção entre fase densa e fase rarefeita.

FIGURA 6:Imagem de simulação por Bernardo Boatini. Sistema em separação de fase no estado estacionário. L=36

Por fim, temos uma imagem do sistema(FIG.7) em ponto crítico $\sigma =0.305$ , onde qualitativamente constata-se que a fase densa e a fase rarefeita não são distinguíveis no sistema e ambas não encontram uma região de estabilidade ao longo da simulação.

FIGURA 7:Imagem de simulação por Bernardo Boatini. Sistema para

\sigma

crítico no estado estacionário. L=36

Para realizar a comparação dos fatores de escala 10 amostras independentes(amostradas a cada $t_{relax}$ ) foram medidas em 4 tamanhos de rede diferentes L = {18, 24, 30, 36}, em cada uma das 10 amostras, 4 subcaixas mediram o número de partículas em diferentes estados estacionários(FIG.8).

FIGURA 8: Compressibilidade Isotérmica(

\mathrm {X} _{L}

) como função de L no sistema de partículas ativas(azul) e para o Ising 2D(laranja)^[1]. Gráfico em escala log-linear com barras de erro iguais ao maior erro padrão extraído dentre os 4 pontos.

É possível ver, que com uma precisão de $R^{2}=0.9805$ a tendência de escala do $\mathrm {X} _{L}$ do sistema de partículas ativas(aproximadamente $\alpha L^{1.4714}$ ) é similar a lei de potência do $\mathrm {X} _{L}$ no Ising 2D.

Os dados nas proximidades do ponto crítico foram extraídos, mas os valores de $U_{4}$ resultantes não foram satisfatórios. Aparentemente a definição de $\langle .\rangle _{L}$ é diferente para as 2 medidas propostas, e por esse motivo a analise de dados foi prejudicada; ou seria necessário uma amostragem com mais de 10 amostra para desprezar o ruído nas medidas, já que em ^[1] o numero de estados estacionários utilizados na estatística é da ordem de $10^{3}$ .

Código

#include <stdio.h>
#include <stdlib.h>
#include <time.h>
#include <math.h>
#include "mc.h"

#define L 30	//L = 18, 24, 30, 36 no papper
#define h 2	//h = 2 no papper
#define c 6	//c = 6 no papper
#define S 10800	//S = (h*L)*(c*L)
#define N 5638	//rho = N/S = 0.522 no papper
int ngb[S][6];
int sb[h*2][L*L];

void sb_creator(void) {

  int i, j, n, aux;
  	//cria subcaixas

  for(n=0; n<h*2; n++){
  	j = 0;
  	aux = n*c*L*L + (c*L/2 - L/2 + j*c*L);
  	if(n<h){
  		for(i=0; i<L*L; i++){
			if(i == (j+1)*L){
				j = j+1;
			}
			aux = n*c*L*L + (c*L/2 - L/2 + j*c*L);
			sb[n][i] = aux + i - j*L;
		}
  	}
  	else{
	aux = h*2 - (n+1);
  		for(i=0; i<L*L; i++){
                        if(i == (j+1)*L){
                        	j = j+1;
                        }
  			if(i-j*L < L/2){
				sb[n][i] = sb[aux][i] + c*L/2;
			}
  			if(i-j*L >= L/2){
				sb[n][i] = sb[aux][i] - c*L/2;
  			}
  		}
  	}
  }

}

void ngb_creator(void) {

  int  i, j, n;

        //cria matriz de vizinhos

  ngb[0][0]= S - c*L + 1;       //superior direita
  ngb[0][1]= 1;                 //direita
  ngb[0][2]= c*L + 1;           //inferior direita
  ngb[0][3]= c*L;               //inferior esquerda
  ngb[0][4]= c*L - 1;           //esquerda
  ngb[0][5]= S - c*L;             //superior esquerda

  for(i=1;i<c*L-1;i++){
        ngb[i][0]= i + S - c*L + 1;
        ngb[i][1]= i + 1;
        ngb[i][2]= i + c*L + 1;
        ngb[i][3]= i + c*L;
        ngb[i][4]= i - 1;
        ngb[i][5]= i + S - c*L;
  }

  ngb[c*L-1][0]= S - c*L;
  ngb[c*L-1][1]= 0;
  ngb[c*L-1][2]= c*L;
  ngb[c*L-1][3]= 2*c*L - 1;
  ngb[c*L-1][4]= c*L - 2;
  ngb[c*L-1][5]= S - 1;

  j = 1;
  n = 2;
  for(i=c*L;i<S-c*L;i++){
        if(i == j*c*L){
                j = j + 1;
                if(j%2 == 0){
                        ngb[i][0]= i - c*L;
                        ngb[i][1]= i + 1;
                        ngb[i][2]= i + c*L;
                        ngb[i][3]= i + 2*c*L - 1;
			ngb[i][4]= i + c*L - 1;
			ngb[i][5]= i - 1;
                }
                else{
                        ngb[i][0]= i - c*L + 1;
                        ngb[i][1]= i + 1;
                        ngb[i][2]= i + c*L + 1;
                        ngb[i][3]= i + c*L;
                        ngb[i][4]= i + c*L - 1;
                        ngb[i][5]= i - c*L;
                }
        }
        else if(i == n*c*L - 1){
                n = n + 1;
                if(n%2 == 0){
                        ngb[i][0]= i - 2*c*L + 1;
                        ngb[i][1]= i - c*L + 1;
                        ngb[i][2]= i + 1;
                        ngb[i][3]= i + c*L;
                        ngb[i][4]= i - 1;
                        ngb[i][5]= i - c*L;
                }
                else{
                        ngb[i][0]= i - c*L;
                        ngb[i][1]= i - c*L + 1;
                        ngb[i][2]= i + c*L;
                        ngb[i][3]= i + c*L - 1;
                        ngb[i][4]= i - 1;
                        ngb[i][5]= i - c*L -1;
                }
        }
        else{
                if(j%2 == 0){
                        ngb[i][0]= i - c*L;
                        ngb[i][1]= i + 1;
                        ngb[i][2]= i + c*L;
                        ngb[i][3]= i + c*L - 1;
                        ngb[i][4]= i - 1;
                        ngb[i][5]= i - c*L - 1;
                }
                else{
                        ngb[i][0]= i - c*L + 1;
                        ngb[i][1]= i + 1;
                        ngb[i][2]= i + c*L + 1;
                        ngb[i][3]= i + c*L;
                        ngb[i][4]= i - 1;
                        ngb[i][5]= i - c*L;
                }
        }
  }

  j = j + 1;

  ngb[S-c*L][0]= (S - c*L) - c*L;
  ngb[S-c*L][1]= (S - c*L) + 1;
  ngb[S-c*L][2]= 0;
  ngb[S-c*L][3]= c*L - 1;
  ngb[S-c*L][4]= (S - c*L) + c*L - 1;
  ngb[S-c*L][5]= (S - c*L) - 1;

  for(i=(S-c*L+1);i<S-1;i++){
        ngb[i][0]= i - c*L;
        ngb[i][1]= i + 1;
        ngb[i][2]= i - (j - 1)*c*L;
        ngb[i][3]= i - (j - 1)*c*L - 1;
        ngb[i][4]= i - 1;
        ngb[i][5]= i - c*L - 1;
  }

  ngb[S-1][0]= (S-1) - c*L;
  ngb[S-1][1]= (S-1) - c*L + 1;
  ngb[S-1][2]= c*L - 1;
  ngb[S-1][3]= c*L - 2;
  ngb[S-1][4]= (S-1) - 1;
  ngb[S-1][5]= (S-1) - c*L - 1; 

return ;
}

int main(void) {

  int i, j, n, nn, r, REP, destiny, MCS_MAX, index, rand, PP[N][2], Space[S], Space_CM[S], d, sum[4];
  double Prob, sigma, randf, x, y, aux, CM;
  REP = 10;
  Prob = 24.0/30.0;
  sigma = 0.305;
  MCS_MAX = 500000;
//-----------------------------------------------INICIA O RNG---------------------------------------------------------
  srand48(time(NULL));
  seed = start_randomic();

//-----------------------------------------------MATRIZ DE VIZINHOS---------------------------------------------------
  ngb_creator();
  sb_creator();

//-----------------------------------------------CRIA O ESTADO INICIAL------------------------------------------------
  for(i=0;i<S;i++){
  	Space[i] = 0;
  	Space_CM[i] = 0;
  }
  for(i=0;i<N;i++){
  	index = S*drand48();
    	if(Space[index]==0){
    		PP[i][0] = index;
		Space[index] = 1;
  	}
	else{
		i = i - 1;
	}
  }
  for(i=0;i<N;i++){
  	index = 6*drand48();
  	PP[i][1] = index;
  }

//-----------------------------------------------DINÂMICA DE MONTE CARLO------------------------------------------------
 for(r=0;r<REP;r++){
 for(n=0;n<MCS_MAX;n++){
  	//printf("\nTempo = %d\n", n+1);
        for(nn=0;nn<N;nn++){
                index = N*drand48();
		rand = round(sigma*ngaussian()); //checar histograma

		// muda polaridade
		if(PP[index][1] + rand > 5){
			PP[index][1] = PP[index][1] + rand - 6;
		}
		else if(PP[index][1] + rand < 0){
			PP[index][1] = PP[index][1] + rand + 6;
		}
		else{
			PP[index][1] = PP[index][1] + rand;
		}
		
		// movimenta
		
		randf = drand48();

		if(Prob >= randf){
			destiny = ngb[PP[index][0]][PP[index][1]];
			if(Space[destiny] == 0){
				Space[PP[index][0]] = 0;
				PP[index][0] = destiny;
				Space[destiny] = 1;
			}
		}
		else{
			rand = 6*drand48();
			destiny = ngb[PP[index][0]][rand];
                        if(Space[destiny] == 0){
                                Space[PP[index][0]] = 0;
                                PP[index][0] = destiny;
                                Space[destiny] = 1;
                        }
		}       
         }
  }
  
  
//-----------------------------------------------CALCULA CENTRO DE MASSA------------------------------------------------
  aux = 0;
  j = 0;
  for(i=0;i<S;i++){
	if(i == (j+1)*c*L){
		j = j + 1;
	}
	if(Space[i] == 1){
  		aux = aux + (i - j*c*L);
	}
  }
  CM = aux/N;
  //printf("Achei o CM=%f\n", CM);
  
//-----------------------------------------------DESLOCA O SISTEMA------------------------------------------------------
  destiny = round(c*L/2);
  index = CM;
  d = destiny - 1 - index;
  //printf("vou deslocar d=%d\n", d);
  for(i=0;i<S;i++){
	if(Space[i]==1){
		destiny = i;
  		for(j=0; j*j<d*d; j++){
			if(d<0){
				destiny = ngb[destiny][4];
			}
			if(d>0){
				destiny = ngb[destiny][1];
			}
  		}	
  		Space_CM[destiny] = 1;
  	}
  }
  
//-----------------------------------------------CONTA PARTICULAS NAS SUBCAIXAS--------------------------------------------

  for(i=0; i<h*2; i++){
  	sum[i] = 0;
  }

  for(i=0; i<h*2; i++){
	for(j=0; j<L*L; j++){
  		sum[i] = sum[i] + Space_CM[sb[i][j]];
	}
  }
  
  printf("%d\t%d\t%d\t%d\n", sum[0], sum[1], sum[2], sum[3]);

/*	
//-----------------------------------------------GERA COORDENADAS PARA PLOT------------------------------------------------

  for(i=0;i<h*L;i++){
  	aux = (1 + pow(-1,i))/4.0;
        for(j=i*c*L;j<(i+1)*c*L;j++){
		y = h*L - i*0.866025403;
		x = 1.0*j - i*c*L + aux;
		if(Space[j]==1){ 
             		printf("%f\t%f\n", x, y);
			}
        }
  }
*/

/*
//-----------------------------------------------CHECA MATRIZ DE PARTICULAS------------------------------------------------
  printf("Site\tPolarity\n");
  for(i=0;i<N;i++){
  	printf("%d\t%d\n", PP[i][0], PP[i][1]);
  }
*/

/*
//-----------------------------------------------CHECA MATRIZ ESPACIAL-----------------------------------------------------
printf("normal\n");
        nn = 0;
        for(i=0;i<h*L;i++){
                for(j=nn*c*L;j<(nn+1)*c*L;j++){
                        if(j == ((nn+1)*c*L)-1){
                                printf(" %d\n", Space[j]);
                        }
                        else{
                                printf(" %d", Space[j]);
                        }
                }
                nn = nn + 1;
        }


//-----------------------------------------------CHECA MATRIZ CENTRALIZADA------------------------------------------------
printf("\ncentralizada\n");
                // checa matriz espacial no CM
        nn = 0;
        for(i=0;i<h*L;i++){
                for(j=nn*c*L;j<(nn+1)*c*L;j++){
                        if(j == ((nn+1)*c*L)-1){
                                printf(" %d\n", Space_CM[j]);
                        }
                        else{
                                printf(" %d", Space_CM[j]);
                        }
                }
                nn = nn + 1;
        }
*/

/*

//-----------------------------------------------CHECA MATRIZ DE VIZINHOS------------------------------------------------ 
  for (i=0; i<S; i++){
  	printf("%d|%d\t%d\t%d\t%d\t%d\t%d\n", i, ngb[i][0], ngb[i][1], ngb[i][2], ngb[i][3], ngb[i][4], ngb[i][5]);  
  }
*/

/*
//-----------------------------------------------CHECA MATRIZ DE SUBCAIXAS-----------------------------------------------
  for(n=0; n<h*2; n++){
  	printf("%d|", n);
  	for(i=0; i<L*L; i++){
        	printf("%d ", sb[n][i]);
  	}
  	printf("\n");
  }
*/
  for(i=0;i<S;i++){
  	Space_CM[i] = 0;
  }
  }
return 0;
}

Referência

↑ ^1,00 ^1,01 ^1,02 ^1,03 ^1,04 ^1,05 ^1,06 ^1,07 ^1,08 ^1,09 ^1,10 ^1,11 ^1,12 ^1,13 ^1,14 https://arxiv.org/pdf/1810.06112.pdf Benjamin Partridge, Chiu Fan Lee, "CRITICAL MOTILITY-INDUCED PHASE SEPARATION BELONGS TO THE ISING UNIVERSALITY CLASS"
↑ Disponível em: https://arxiv.org/pdf/0709.1325.pdf Alexandre Lefèvre, Giulio Biroli, "DYNAMICS OF INTERACTING PARTICLE SYSTEMS: STOCHASTIC PROCESS AND FIELD THEORY"
↑ ^3,0 ^3,1 HATZIKIROU, Haralambos et al. Extracting cellular automaton rules from physical Langevin equation models for single and collective cell migration. Journal of Mathematical Biology, [s. l.], 2017.

[BENJAMIN-1] 1,00 ^1,01 ^1,02 ^1,03 ^1,04 ^1,05 ^1,06 ^1,07 ^1,08 ^1,09 ^1,10 ^1,11 ^1,12 ^1,13 ^1,14 https://arxiv.org/pdf/1810.06112.pdf Benjamin Partridge, Chiu Fan Lee, "CRITICAL MOTILITY-INDUCED PHASE SEPARATION BELONGS TO THE ISING UNIVERSALITY CLASS"

[ALEXANDRE-2] Disponível em: https://arxiv.org/pdf/0709.1325.pdf Alexandre Lefèvre, Giulio Biroli, "DYNAMICS OF INTERACTING PARTICLE SYSTEMS: STOCHASTIC PROCESS AND FIELD THEORY"

[CELULAR-3] 3,0 ^3,1 HATZIKIROU, Haralambos et al. Extracting cellular automaton rules from physical Langevin equation models for single and collective cell migration. Journal of Mathematical Biology, [s. l.], 2017.

[1]

[2]

[3]

Motility-Induced Phase Separation(MIPS)

Índice

Introdução