---------------------------------------
TAREA COLECTIVA:
Encontrar la dependencia de la fraccion de volumen
de polimeros critica (etap_critica) con 
el cociente de radios polimero-coloide (RATIO)
Bibliografia:
http://www.uam.es/personal_pdi/ciencias/rdelgado/FELSC.html
-Apuntes del curso 09/10,
-Apuntes Coloides paginas 7-9.


UNIDADES
Longitud: diametro del coloide.
energia: k_B T (el potencial no depende de la temperatura usando esta unidad)

PARAMETROS: NP,Lx,Ly,Lz,ETAP,RATIO
Fraccion de ocupacion de coloides-> ETAC=(pi/6.)*Nparticulas/(lx*ly*lz)
Fraccion de ocupacion de coloides-> ETAP
Cociente RADIO_POLIMERO/RADIO_COLOIDE-> RATIO 

---------------------------------------------
POTENCIAL DE DEPLECION ASAKURA-OOSAWA: 
ener(dist)
dist=distancia entre dos coloides
defino: 
ratm1=ratio+1
ddd=dist/ratm1
rrr=(ratm1/ratio)**3
...
ener(dist)=INFINITO                          SI DIST < 1  (SOLAPAMIENTO ESFERAS DURAS)
ener(dist)=etap*rrr*(1.5*ddd-0.5*ddd**3-1.0) SI 1< DIST < RATM1
ener(dist)=0                                 SI DIST>RATM1
...

Nota: Para implementar solapamientos en nucleo duro
puede usarse una variable logica "accept" (ver energy.f en mc_ao.tar -pag.web-)
---------------------------------------------

OBJETIVO:

Encontrar la dependencia de la fraccion de volumen
de polimeros critica (etap_critica) con 
el cociente de radios RATIO=R_POLYMER/R_COLOIDE

METODO
Sabemos que la concentracion critica de coloides es cercana a ETAC=0.15.
Asi que trabajamos con esta fraccion de volumen de coloides fija 
y con un valor fijo de RATIO.
Variamos ETAP hasta encontrar la frontera de 
coexistencia de fases  (liquido_coloide  + gas_coloide).
En la frontera hay una transicion de una fase homogenea a "dos fases".


CADA ALUMNO TIENE ASIGNADO UN VALOR DE "RATIO"
entre RATIO=0.1 y RATIO=0.8 (se asignara en clase)
Por supuesto se pueden hacer varios "ratios", si uno quiere.
---------------------------------------------

PUESTA A PUNTO de Tu Codigo

Si vais a usar vuestro codigo, verificad primero que la energía por particula
que obteneis el la misma que esta:

ETAC         ETAP     RATIO      ENERGY      LBOX
0.135        0.2       0.2       -0.567       14.9         
0.2618       0.4       0.2       -6.07        10.0

Notas: 
      Podeis usar fcc.f (www.uam.es/rafael.delgado/docencia.html...)
      para generar la configuracion inicial al valor de ETAC deseado
      Podeis cojer otro Lbox (lado caja), eso puede dar pequeñas desviaciones
      en la energia.
      La energía por particula es la mitad de la energia de interaccion 
      de una particula con el resto, pues las interacciones son a pares
      (es decir, una mitad para la particula y la otra para el resto).
      El valor de ENERGY en la tabla de arriba es DOS VECES LA ENERGIA POR PARTICULA.

---------------------------------------

METODO OPERACIONAL

Nota: Pondremos siempre todas las particulas en la caja 
(es decir, condiciones periodicas en las posiciones)	

1) Generar configuracion inicial usando fcc.f al valor de ETAC=0.15
(el fichero output es fcc.rst, es leido por mc.f en la subroutina init.f)
conviene no elegir LBOX muy pequeño (mejor contra mas grande).
Por ejemplo LBOX>10, para evitar en lo posible efectos de tamaño finito

2) Elegid un valor de ETAP 
(nota: normalmente para ETAP<0.2 la fase es homogenea;
para ETAP mayores se encuentra la transicion)

---------------------------------------

PRIMER PASO SIMULACION
   
Corred el codigo a partir del fichero con posiciones fcc hasta
alcanzar el equilibrio. Para determinar si estais en el equilibrio,
podeis pintar la energia. En el equilibrio fluctuara en torno 
a una media bien definida.


Guardad
a)   e=E/N, energia por particula 
b)   func. distr. radial  (la subrutina va adjunta en mc_ao.tar descargable desde 
			   mi pag. web)
c)   fichero de posiciones final,
     que sera usado posteriormente como "restart file", 
     en el codigo mc.f se llama "mc_save.rst".

---------------------------------------

SEGUNDO PASO SIMULACION

   -Cambiad el tamaño de la caja en el restart file
    Aumentad LX, Por ejemplo LX->LX+20.0

    -Guardad
   a)   Varianza de las posiciones en direccion x.
        Var_x= (1/N) sum_i (xi-xmedio)**2		
        (podeis hacer xmedio=0, si la caja esta centrada en cero)	

    Si el sistema se separa en dos fases, aparece una interfase
    liquido-gas, y Var_x no aumentara mucho.
    Pero en una fase homogenea la varianza subira rápidamente 
    al ocupar el volumen extra que le hemos dado aumentando LX.	
     
   b) Perfil de densidad (opcional) 
      habria que usar vuestra subrutina histograma.
      Es interesante para calcular la densidad de equilibrio de cada fase.

---------------------------------------

TERCER PASO
Decidid si el sistema es homogeneo o esta dentro de la curva de coexistencia
(en este caso la fase liquida no se expande por toda la caja).
-Para ello podeis usar el ejecutable "punto" (en mi pag. web) para ver las particulas
 y/o mirar la señal de Var_x

Si es homogeneo, REPETIR PRIMER Y SEGUNDO PASO, aumentando ETAP,
si es bifasico, REPETIR, con el promedio entre el valor ETAP homogeneo
mas alto que tengais y el valor ETAP bifasico mas pequeño.
De este modo podeis estimar ETAP_critico, al menos con precision 
ERROR_ETAP = +-0.05.
-Estimo que seran necesarias unas 5 simulaciones, mas o menos.-