Experimento en favor de la teoría de Jeremy England (III)

Encendemos un termostato sobre un recipiente de volumen constante con un gas monoatómico en su interior siguiendo un potencial Lennard-Jones.

Dicho potencial tienen una componente positiva que resulta en una fuerza repulsiva actuando a pequeñas distancias entre partículas (simulando el resultado de la sobreposición de los orbitales electrónicos, conocido como la repulsión de Pauli). De esta manera, las partículas no puede solaparse a partir de cierta distancia mínima. Por otra parte, hay una componente negativa del potencial que simboliza una fuerza atractiva que actúa a mayores distancias entre partículas (simulando la fuerza de Van Der Waals, o fuerza de dispersión).

Comportamiento normal del sistema en estas condiciones.

Si comenzamos con un sistema de N partículas monoatómicas inicialmente en reposo, al activar el termostato a una determinada temperatura kT, el sistema comienza a adquirir energía cinética y las partículas van adquiriendo velocidad. El potencial L-J descrito decide más tarde cómo se producen las fuerzas interacción entre partículas (evitando solapamiento, y favoreciendo la atracción mutua entre pares en relación inversa a su distancia).

Si el termostato se mantiene fijo en la misma temperatura kT, pasado un cierto número de segundos el sistema alcanza un equilibrio térmico caótico de máxima entropía: en este punto no se observan ya grandes variaciones en las medidas energéticas y de velocidad.

Esta sucesión de acontecimientos es la que se observa siempre que se repite el experimento no importa el número de veces que se realice, ni el tiempo que se espere. El hecho de que un fenómeno complejo ocurra en estas condiciones, como por ejemplo que N / 2 partículas se agrupen muy cerca unas de otras es casi imposible dada las pocas configuraciones que muestran este estado en comparación con los billones de configuraciones que presentan un estado caótico como el observado.



El poder de la computación evolutiva.

No voy a repetir el modo en que funciona la computación evolutiva, porque ya lo he hecho en varias entradas anteriores (por ejemplo en esta, y también aquí).

Baste comentar que mediante computación evolutiva podemos simular cientos de miles de sistemas físicos como el anteriormente descrito, y que además podemos buscar de un modo muy eficiente cualquier comportamiento deseado: por ejemplo, podemos pretender encontrar aquellos sistemas que, partiendo de un mismo estado inicial, evolucionen hacia estados que tengan mucha energía cinética, o hacia aquellos que tengan poca separación entre partículas, etc.

Pues bien, en un sistema como el descrito hasta ahora, la probabilidad de que el camino seguido por cualquiera de los miles de sistemas termine en un estado con cierto orden (con partículas muy cercanas ocupando poco espacio del volumen disponible, por ejemplo), es tan pequeña que no se llega nunca a presentar. Por mucho que dejemos pasar el tiempo de computación, y por mucho que seleccionemos en cada iteración los sistemas que más cercanía media presentan entre partículas, el resultado final apenas se diferencia de un estado caótico.

La evolución de este tipo de sistemas está determinado, en la práctica, a presentar sistemas caóticos.

¿Cómo conseguir mejorar la probabilidad de encontrar orden en el sistema?

Pues hasta hace poco no se tenía mucha idea, pero desde hace unos años para acá, ciertos físicos de renombre están intentando arrojar luz sobre el asunto. Uno de ellos es Jeremy England, físico del MIT en Estados Unidos.

Su propuesta se resumen en la siguiente fórmula (que es intimidadora en un primer vistazo, pero una vez se entiende cada componente de la misma no lo es tanto):

Lo que este autor viene a decir, resumiendo mucho, es que un sistema como el que hemos descrito en el apartado anterior no puede evolucionar hacia estados ordenados (disminuyendo la entropía del sistema) a menos que se produzca una buena absorción y consumo constante de nueva energía externa (convirtiéndose en un sistema abierto). Es decir, a menos que se produzca una generación de calor extra gracias a este mayor consumo energético.

La probabilidad de que un sistema ordenado acontezca va a estar pues directamente relacionada con la cantidad de energía externa consumida por el sistema. Complejidad y consumo energético van ligados de la mano en estos sistemas abiertos, y aquellos sistemas que durante su recorrido en el tiempo muestran una mayor eficiencia recibiendo esta energía externa (aunque sea de modo fortuito), serán los que presentarán con mayor probabilidad algo de orden (incremento local de entropía negativo).

Simulación de este sistema abierto.

Para experimentar con esta propuesta, vamos simplemente a añadir una nueva fuente externa de energía, que simularemos a través de colisiones arbitrarias con algunas de las partículas. En concreto, vamos a simular una lluvia de "fotones", que cada cierto tiempo, con una cierta probabilidad, colisionará con alguna de nuestra N partículas, simulándose en este caso la absorción del "fotón" con el consiguiente calor disipado y el cambio en la dirección de la partícula.

Es decir, que hemos simulado una lluvia aleatoria de "fotones" entrando y saliendo del sistema, los cuales pueden fortuitamente a veces colisionar con una de nuestras partículas, produciéndose su absorción por nuestro átomo y la generación de calor implícita a este proceso, además del cambio en el momento de la partícula (cambio en su vector de velocidad).

Comportamiento normal del sistema en estas nuevas condiciones.

Si dejamos el nuevo sistema abierto proceder de modo espontáneo durante un tiempo, parece que poco ha cambiado. El sistema comienza a adquirir velocidad por el termostato, y finalmente se llega a un estado caótico de equilibrio muy similar al del caso anterior con el sistema cerrado (sin lluvia de "fotones").

Sin embargo, esto es sólo una falsa impresión. El nuevo sistema es radicalmente distinto en un importante detalle: ahora la probabilidad de alcanzar un estado ordenado ha aumentado mucho, no siendo ya tal orden un caso fortuito de unas pocas configuraciones posibles entre billones, sino que la probabilidad se ve mejorada por la tercera y cuarta componente de la fórmula arriba vista de Jeremy. La probabilidad de ver orden en un sistema abierto aumenta para aquellos caminos que más y mejor consumen (y disipan en forma de calor) esta nueva energía externa (aunque sea de un modo fortuito).

Si a este nuevo sistema abierto le aplicamos una selección por computación evolutiva, ahora la probabilidad de ciertas configuraciones es lo suficientemente alta como para que lleguen a acontecer gradualmente al pasar las iteraciones (el tiempo).

Más complejidad requiere más consumo, y más consumo requiere más complejidad.

La lluvia de "fotones" no sólo pone al descubierto la teoría principal de Jeremy England de que en sistemas abiertos a una fuente de energía externa ciertas configuraciones complejas pueden ser lo suficientemente probables como para que gradualmente se acumulen y crezcan en complejidad, sino que también pone de evidencia el hecho de que cuanto más energía absorbe el sistema, más complejidad puede alcanzar y mantener, y el hecho inverso; que cuanto más complejo es un sistema que se vea aparecer, más energía debe haber consumido en el tiempo transcurrido (alta eficiencia).

Metodología del experimento.

No voy a entrar en mucho detalle porque, como digo, ya he explicado bastante a fondo el asunto en anteriores entradas, pero sí quiero hacer notar, que este nuevo artículo viene al cuento  porque he modificado todo el código fuente del programa para adaptarlo a entornos 3D (hasta ahora, estas pruebas las hacía en 2D).

He desarrollado un código totalmente nuevo en un Applet de Java,  partiendo en parte del anterior en dos dimensiones, y también de parte de un desarrollo de código libre ya existente y desarrollado por Cheng Zhang. El código no está excesivamente comentado, por lo que cualquier aportación que comente y limpie un poco el código será bienvenida ;).

Cosas por hacer.

Pues hay infinidad de mejoras y ampliaciones posibles, y por desgracia no tengo tiempo libre para hacer todo lo que quisiera, pero una cosa que sí deseo programar, es un cambio que permita observar otro tipo de patrones complejos como vórtices, torbellinos, y cosas así. La idea en principio es clara: colocar algunos átomos fijos e inamovibles en algún plano (u otra configuración), y ver cómo responde el sistema con y sin lluvia de "fotones" (sistema abierto o cerrado).

Edito. 24/11/2015

He añadido, como ha aconsejado Sergio en un comentario, dos tipos de partículas distintas (rojas y verdes), las cuales modifican el potencial L-J de base del siguiente modo: las partículas iguales se repelen un poco entre sí (en valor inverso a su distancia), mientras que las partículas distintas se atraen (también en relación inversa con su distancia). El resultado es bastante interesante:

Ejemplo del nuevo sistema de dos partículas alcanzando el equilibrio (caos)

Cúmulo de agrupaciones que aparece cuando se buscan aquellos sistemas que más energía absorben (y disipan en forma de calor)

Este experimento muestra de nuevo la estrecha relación que existen en los sistemas naturales entre energía y complejidad estructural (o espacial). En concordancia con la teoría de Jeremy Enalgand, se puede observar que un gran aumento en la irreversibilidad del sistema (un gran orden como el mostrado en esta última imagen) debe ir siempre acompañado del consumo de una cantidad de energía muy superior a la media de un sistema caótico (y por lo tanto bastante más fácilmente reversible).

Anexo.

Os dejo a continuación a modo de anexo todo el código fuente desarrollado hasta ahora en esta versión 3D, por si alguien desea continuar por este camino:

1) LJ3MDApp.java.

package md_3d;

/* Simple molecular dynamics simulation of a Lennard-Jones fluid

Copyright 2011-2014 Cheng Zhang
 
Modified by Samuel Graván 2015

This program is free software: you can redistribute it and/or modify
it under the terms of the GNU General Public License as published by
the Free Software Foundation, either version 2 of the License, or
(at your option) any later version.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
GNU General Public License for more details.

A copy of the GNU General Public License can be found in
<http://www.gnu.org/licenses/>.
*/


import static java.lang.Math.PI;
import static java.lang.Math.pow;
import static java.lang.Math.sqrt;

import java.applet.Applet;
import java.awt.BorderLayout;
import java.awt.Color;
import java.awt.Component;
import java.awt.Container;
import java.awt.Dimension;
import java.awt.Font;
import java.awt.Graphics;
import java.awt.GridLayout;
import java.awt.Image;
import java.awt.event.ActionEvent;
import java.awt.event.ActionListener;
import java.awt.event.MouseEvent;
import java.awt.event.MouseListener;
import java.awt.event.MouseMotionListener;
import java.awt.event.MouseWheelEvent;
import java.awt.event.MouseWheelListener;
import java.awt.image.IndexColorModel;
import java.awt.image.MemoryImageSource;
import java.text.DecimalFormat;
import java.util.Random;

import javax.swing.JApplet;
import javax.swing.JButton;
import javax.swing.JCheckBox;
import javax.swing.JLabel;
import javax.swing.JPanel;
import javax.swing.JTextField;
import javax.swing.JToggleButton;
import javax.swing.Timer;

public class LJ3MDApp extends JApplet implements ActionListener {
 double rhoInit = 0.005;
 LJ3MD md = new LJ3MD(rhoInit);
 int delay = 100;
 int speed = 10; // md steps per frame
 Timer timer;

 XYZCanvas canvas; // 3D drawing here
 JPanel cpnl, spnl;
 JTextField tNum = new JTextField("" + md.N);
 JTextField tTemp = new JTextField("15.0");
 JTextField tRho = new JTextField("" + rhoInit);
 JTextField tSpeed = new JTextField("" + speed);
 JButton bReset = new JButton("Reset");
 JButton bRetime = new JButton("Reset time");
 JToggleButton bStart = new JToggleButton("Start");
 JToggleButton bEnvolve = new JToggleButton("Evolution");
 JCheckBox bTstat = new JCheckBox("Thermostat", true);
 JLabel lStatus = new JLabel("Status");
 JTextField tAvK = new JTextField("0");
 JTextField tAvU = new JTextField("0");
 JTextField tAvp = new JTextField("0");
 public static final long serialVersionUID = 1L;

 public void init() {
  tNum.setHorizontalAlignment(JTextField.CENTER);
  tTemp.setHorizontalAlignment(JTextField.CENTER);
  tRho.setHorizontalAlignment(JTextField.CENTER);
  tSpeed.setHorizontalAlignment(JTextField.CENTER);

  tAvK.setHorizontalAlignment(JTextField.RIGHT);
  tAvU.setHorizontalAlignment(JTextField.RIGHT);
  tAvp.setHorizontalAlignment(JTextField.RIGHT);

  Container box = getContentPane();
  box.setLayout(new BorderLayout());

  cpnl = new JPanel(); // create a panel for controls
  cpnl.setLayout(new GridLayout(19, 2));
  box.add(cpnl, BorderLayout.EAST);

  // add controls
  cpnl.add(bStart);
  bStart.addActionListener(this);

  cpnl.add(bEnvolve);
  bEnvolve.addActionListener(this);

  cpnl.add(bReset);
  bReset.addActionListener(this);

  cpnl.add(new JLabel(" N:"));
  tNum.addActionListener(this);
  cpnl.add(tNum);

  cpnl.add(new JLabel(" Temperature:"));
  tTemp.addActionListener(this);
  cpnl.add(tTemp);

  cpnl.add(new JLabel(" Density (\u03c1):"));
  tRho.addActionListener(this);
  cpnl.add(tRho);

  cpnl.add(new JLabel(" Steps/frame:"));
  tSpeed.addActionListener(this);
  cpnl.add(tSpeed);

  cpnl.add(bTstat);
  bTstat.addActionListener(this);

  cpnl.add(new JLabel(" < K/N > :"));
  tAvK.setEditable(false);
  cpnl.add(tAvK);

  cpnl.add(new JLabel(" < U/N > :"));
  tAvU.setEditable(false);
  cpnl.add(tAvU);

  cpnl.add(new JLabel(" < pressure > :"));
  tAvp.setEditable(false);
  cpnl.add(tAvp);

  cpnl.add(bRetime);
  bRetime.addActionListener(this);

  spnl = new JPanel(); // create a panel for status
  lStatus.setPreferredSize(new Dimension(1200, 20));
  box.add(spnl, BorderLayout.SOUTH);
  lStatus.setFont(new Font("Courier", 0, 12));
  spnl.add(lStatus);

  canvas = new XYZCanvas();
  box.add(canvas, BorderLayout.CENTER);

  timer = new Timer(delay, this);
  timer.start();
  timer.stop();
 }

 public void start() {
 }

 public void update(Graphics g) {
  paint(g);
 }

 DecimalFormat df = new DecimalFormat("###0.000");

 // int frame = 0;
 public void paint(Graphics g) {
  // System.out.println("frame: " + (frame++));
  lStatus.setText("t = " + df.format(md.dt * md.step) + ", " + "N = " + md.N + ", " + "E/N = "
    + df.format(md.E / md.N) + ", " + "U/N = " + df.format(md.U / md.N) + ", " + "K/N = "
    + df.format(md.K / md.N) + ", " + "p = " + df.format(md.p) + "; "
    + "avK = " + df.format(md.avK.getAve())
    + ", avU = " + df.format(md.avU.getAve())
    + ", avE = " + df.format(md.avK.getAve() + md.avU.getAve()));
  tAvK.setText(df.format(md.avK.getAve() / md.N) + "  ");
  tAvU.setText(df.format(md.avU.getAve() / md.N) + "  ");
  tAvp.setText(df.format(md.avp.getAve()) + "  ");
  canvas.refresh(md.getXWrap(), md.N, true, false);
  cpnl.repaint();
  spnl.repaint();
  lStatus.repaint();
 }

 public void actionPerformed(ActionEvent e) {
  Object src = e.getSource();
  if (src == timer) {
   for (int i = 0; i < speed; i++)
   {
    if(md.step % 100 == 0 && md.step > 0)
    {
     md.vvEvolution();
     md.step++;
    }
    else
    {
     md.vv();
    }
   }
   repaint();
   return;
  }

  boolean adjCanvasScale = false;

  if (src == bTstat)
   md.thermostat = !md.thermostat;

  if (src == tTemp || src == bReset) {
   double kT = Double.parseDouble(tTemp.getText().trim());
   if (kT < 1e-8) {
    kT = 1e-8;
    tTemp.setText("  " + kT);
   }
   md.kT = kT;
   md.clearData();
  }

  if (src == tRho || src == bReset) {
   double rho = Double.parseDouble(tRho.getText().trim());
   if (rho < 1e-3) {
    rho = 1e-3;
    tRho.setText("   " + rho);
   }
   if (rho > 1.2) {
    rho = 1.2;
    tRho.setText("   " + rho);
   }
   md.setDensity(rho);
   md.clearData();
   adjCanvasScale = true;
  }

  if (src == tSpeed || src == bReset) {
   speed = Integer.parseInt(tSpeed.getText().trim());
   if (speed < 1) {
    speed = 1;
    tSpeed.setText("   " + speed);
   }
  }

  if (src == bRetime)
   md.clearData();

  if (src == bStart) {
   boolean on = bStart.isSelected();
   if (on) {
    timer.restart();
    bStart.setText("Pause");
   } else {
    timer.stop();
    bStart.setText("Resume");
   }
  }

  if (src == bEnvolve) {
   LJ3MD mdI = new LJ3MD(0.005f);

   mdI.init(0.005f);
   for (int k = 0; k <= 100; k++) {
    mdI.vv();
   }

   EvolutionLibrary evolution = new EvolutionLibrary(250, 1250, 100);
   Poblacion resultados = evolution.buscarSelection(mdI, canvas, this, this.getGraphics());

   md.clearData();
   
   mdI = new LJ3MD(0.005f);
   mdI = resultados.getIndividuos().get(0).getLj();

   md.init(mdI.getDen(), mdI.getN(), mdI.getV(), mdI.getF(), mdI.getX(), mdI.getStep(), mdI.getRho(),
     mdI.getDt(), mdI.getVol(), mdI.getRc(), mdI.getkT(), mdI.getXWrap(), mdI.getK(), mdI.getU(),
     mdI.getUs(), mdI.getE(), mdI.getVir(), mdI.getP(), mdI.getUshift(), mdI.getPtail(), mdI.getUtail(),
     mdI.getAvU(), mdI.getAvK(), mdI.getAvp(), mdI.getXij());

   bEnvolve.setText("Done!!");
  }

  if (src == tNum) {
   int n = Integer.parseInt(tNum.getText().trim());
   if (n < 2) {
    n = 2;
    tNum.setText(" " + n);
   }
   md.N = n;
   md.init(md.rho);
   adjCanvasScale = true;
  }

  if (src == bReset) {
   if (timer.isRunning())
    timer.stop();
   bStart.setSelected(false);
   bStart.setText("Start");
   md.init(md.rho);
  }

  canvas.refresh(md.getXWrap(), md.N, true, adjCanvasScale);

  repaint();
 }
}

class LJ3MD {
 static final int D = 3;
 public int N = 64; // number of particles 8^3/2
 private int DOF = N * D - D * (D + 1) / 2;
 double den;
 double rho = 0.256;
 double dt = 0.002; // time step for integrating Newton's equation
 public double L = 10.0; // side of the box
 public double Vol = 1000.0;
 double rc = 8d;
 public double kT = 15.0; // temperature times Boltzmann constant
 double x[][]; // position, from 0 to 1
 double v[][], f[][]; // velocity and force
 double xwrap[][]; // wrapped coordinates
 double K, U, Us, E; // kinetic, potential, and total energy
 double Vir, p; // virial and pressure
 boolean thermostat = true; // use thermostat
 double ushift;
 double Utail; // potential energy tail correction
 double ptail; // pressure tail correction
 int step; // simulation steps
 Ave avU = new Ave(), avK = new Ave(), avp = new Ave();

 /** constructor */
 LJ3MD(double den) {
  init(den);
 }

 /** initialize system */
 public void init(double den) {
  int i, j, k, id;

  this.den = den;
  DOF = N * D - D * (D + 1) / 2;
  setDensity(den);
  step = 0;
  x = new double[N][D];
  v = new double[N][D];
  f = new double[N][D];
  xwrap = new double[N][D];
  for (id = 0; id < N; id++)
   for (j = 0; j < D; j++)
    v[id][j] = 0.0;

  int N1 = (int) (pow(2 * N, 1.0 / D) + .999999); // # of particles per
              // side (cubic lattic)
  double a = 1. / N1;
  for (id = 0, i = 0; i < N1 && id < N; i++) // place particles on a cubic
             // lattice
   for (j = 0; j < N1 && id < N; j++)
    for (k = 0; k < N1 && id < N; k++) {
     if ((i + j + k) % 2 != 0)
      continue;
     x[id][0] = a * (i + .5);
     x[id][1] = a * (j + .5);
     x[id][2] = a * (k + .5);
     id++;
    }
 }

 /** clonate a system */
 public void init(double den, int N, double v[][], double f[][], double x[][], int step, double rho, double dt,
   double Vol, double rc, double kT, double xwrap[][], double K, double U, double Us, double E, double Vir,
   double p, double ushift, double ptail, double Utail, Ave avU, Ave avK, Ave avp, double xij[]) {
  this.rho = rho;
  this.den = den;
  init(den);
  this.N = N;
  this.v = new double[N][D];
  for (int k = 0; k < N; k++) {
   for (int j = 0; j < D; j++)
   {
    this.v[k][j] = v[k][j];
   }
  }
  this.f = new double[N][D];
  for (int k = 0; k < N; k++) {
   for (int j = 0; j < D; j++)
   {
    this.f[k][j] = f[k][j];
   }
  }
  this.x = new double[N][D];
  for (int k = 0; k < N; k++) {
   for (int j = 0; j < D; j++)
   {
    this.x[k][j] = x[k][j];
   }
  }
  this.step = step;  
  this.dt = dt;
  this.Vol = Vol;
  this.rc = rc;
  this.kT = kT;
  this.xwrap = new double[N][D];
  for (int k = 0; k < N; k++) {
   for (int j = 0; j < D; j++)
   {
    this.xwrap[k][j] = xwrap[k][j];
   }
  }
  this.K = K;
  this.U = U;
  this.Us = Us;
  this.E = E;
  this.Vir = Vir;
  this.p = p;
  this.ushift = ushift;
  this.ptail = ptail;
  this.Utail = Utail;
  
  this.avK = new Ave();
  this.avK.cnt = avK.cnt;
  this.avK.x2sm = avK.x2sm;
  this.avK.xsm = avK.xsm;
  
  this.avU = new Ave();
  this.avU.cnt = avU.cnt;
  this.avU.x2sm = avU.x2sm;
  this.avU.xsm = avU.xsm;
  
  this.avp = new Ave();
  this.avp.cnt = avp.cnt;
  this.avp.x2sm = avp.x2sm;
  this.avp.xsm = avp.xsm;

  this.xij = new double[D];
  for (int k = 0; k < D; k++) {
   this.xij[k] = xij[k];
  }
 }

 /** set density and recalculate tail corrections */
 void setDensity(double den) {
  rho = den;
  Vol = N / rho;
  L = pow(Vol, 1.0 / D);

  /* compute shifts and tail corrections */
  //if ((rc = 2.5) > .5 * L)
  // rc = .5 * L;
  double irc = 1. / rc, irc3 = irc * irc * irc, irc6 = irc3 * irc3;
  ushift = 4.0 * irc6 * (irc6 - 1);
  Utail = 8.0 * PI / 9 * rho * N * (irc6 - 3) * irc3;
  ptail = 16.0 * PI / 9 * rho * rho * (irc6 - 1.5) * irc3;

  clearData();
 }

 void clearData() {
  step = 0;
  avU.clear();
  avK.clear();
  avp.clear();
 }

 public double getMeanEnergy() {
  return E / step;
 }

 public double getMeanFreedom() {
  double freedom = 0d;
  double xijF[] = new double[D];

  for (int i = 0; i < N - 1; i++) {
   for (int j = i + 1; j < N; j++) {
    
    double dx = pbc(x[i][0] - x[j][0]);
    double dy = pbc(x[i][1] - x[j][1]);
    double dz = pbc(x[i][2] - x[j][2]);

    freedom += dx * dx + dy * dy + dz * dz;
   }
  }

  return freedom;
 }

 /** integrate Newton's equation by one step, velocity Verlet */
 void vv() {
  double dth = dt * .5, dtl = dt / L;

  for (int id = 0; id < N; id++) // velocity-verlet, first half
   for (int d = 0; d < D; d++) {
    v[id][d] += f[id][d] * dth;
    x[id][d] += v[id][d] * dtl;
   }

  force(); // compute force

  for (int id = 0; id < N; id++) // velocity-verlet, second half
   for (int d = 0; d < D; d++)
    v[id][d] += f[id][d] * dth;

  rmcom(); // rm angular momentum however isn't good idea due to
     // pbc and minimal image convension
  K = vrescale(0.02); // compute kinetic energy and adjust velocity
  E = K + Us;
  step++;
  avU.add(U);
  avK.add(K);
  avp.add(p);
 }

 /** Simulamos colisiones aleatorias entre una partícula del sistema y
  *  una partícula externa que introduce energía en nuestro sistema al ser absorbida,
  *  con el consiguiente cambio en el vector velocidad de nuestra partícual.
  *  La partícual externa desaparece (es absorbida) tras la colisión.  */
 void vvEvolution() {
  for (int id = 0; id < N; id++)
  {
   Random r = new Random();
   double pr = r.nextInt(100);

   if (pr < 10) {
    v[id][0] *= -1;
    v[id][1] *= -1;
    v[id][2] *= -1;
   }
  }
 }

 /**
  * calculate potential energy and force half-box cutoff, minimal distance
  * image
  */
 double xij[] = new double[D];

 void force() {
  int i, j, d, prcnt = 0;
  double invr2, invr6, r2, fscal, tmp;

  // clear the force and energy
  U = 0.0;
  Vir = 0.0;
  for (i = 0; i < N; i++)
   for (d = 0; d < D; d++)
    f[i][d] = 0.0;

  // loop over pairs of particles and compute energy
  for (i = 0; i < N - 1; i++) {
   for (j = i + 1; j < N; j++) {
    rv3Diff(xij, x[i], x[j]);
    vpbc(xij);
    r2 = vsqr(xij);
    if (r2 > rc * rc)
     continue;
    invr2 = 1.0 / r2;
    invr6 = invr2 * invr2 * invr2;
    fscal = 48.0 * invr6 * (invr6 - 2.5);
    Vir += fscal;
    fscal *= invr2;
    for (d = 0; d < D; d++) {
     tmp = xij[d] * fscal;
     f[i][d] += tmp;
     f[j][d] -= tmp;
    }
    U += 4.0 * invr6 * (invr6 - 1.0);
    prcnt++;
   }
  }
  Us = U - prcnt * ushift;
  U += Utail;
  p = rho * kT + Vir / (3 * Vol) + ptail;
 }

 Random rng = new Random();

 /** velocity rescaling thermostat */
 double vrescale(double dt) {
  double ek1 = getEKin();
  if (!thermostat)
   return ek1;
  double ekav = .5f * kT * DOF;
  double amp = 2 * sqrt(ek1 * ekav * dt / DOF);
  double ek2 = ek1 + (ekav - ek1) * dt + amp * rng.nextGaussian();
  if (ek2 < 0)
   ek2 = ek1;
  if (ek1 > 0) {
   double s = sqrt(ek2 / ek1);
   for (int i = 0; i < N; i++)
    for (int d = 0; d < D; d++)
     v[i][d] *= s;
  }
  return ek2;
 }

 /** remove motion of the center of mass */
 void rmcom() {
  for (int d = 0; d < D; d++) { // remove the center of mass motion
   double vc = 0.0;
   for (int id = 0; id < N; id++)
    vc += v[id][d];
   vc /= N;
   for (int id = 0; id < N; id++)
    v[id][d] -= vc;
  }
 }

 private void rv3Diff(double c[], double a[], double b[]) {
  c[0] = a[0] - b[0];
  c[1] = a[1] - b[1];
  c[2] = a[2] - b[2];
 }

 private double rv3Dot(double a[], double b[]) {
  return a[0] * b[0] + a[1] * b[1] + a[2] * b[2];
 }

 /** return the kinetic energy */
 double getEKin() {
  double ek = 0.0;
  for (int i = 0; i < N; i++)
   ek += rv3Dot(v[i], v[i]);
  return .5 * ek;
 }
 
 /** image with the minimal distance */
 double pbc(double a) {
  return L * (a - ((int) (a + 1000.5) - 1000));
 }

 void vpbc(double v[]) {
  v[0] = pbc(v[0]);
  v[1] = pbc(v[1]);
  v[2] = pbc(v[2]);
 }

 double vsqr(double v[]) {
  return v[0] * v[0] + v[1] * v[1] + v[2] * v[2];
 }

 /** wrap coordinates back into box */
 public double[][] getXWrap() {
  for (int i = 0; i < N; i++)
   for (int d = 0; d < D; d++) {
    double xx = x[i][d] + 1000.0;
    xwrap[i][d] = (xx - (int) xx) * L;
   }
  return xwrap;
 }

 public double getK() {
  return K;
 }

 public void setK(double k) {
  K = k;
 }

 public double getU() {
  return U;
 }

 public void setU(double u) {
  U = u;
 }

 public double getE() {
  return E;
 }

 public void setE(double e) {
  E = e;
 }

 public Ave getAvU() {
  return avU;
 }

 public void setAvU(Ave avU) {
  this.avU = avU;
 }

 public Ave getAvK() {
  return avK;
 }

 public void setAvK(Ave avK) {
  this.avK = avK;
 }

 public int getN() {
  return N;
 }

 public void setN(int n) {
  N = n;
 }

 public double getRho() {
  return rho;
 }

 public void setRho(double rho) {
  this.rho = rho;
 }

 public double getDt() {
  return dt;
 }

 public void setDt(double dt) {
  this.dt = dt;
 }

 public double getVol() {
  return Vol;
 }

 public void setVol(double vol) {
  Vol = vol;
 }

 public double getRc() {
  return rc;
 }

 public void setRc(double rc) {
  this.rc = rc;
 }

 public double getkT() {
  return kT;
 }

 public void setkT(double kT) {
  this.kT = kT;
 }

 public double[][] getX() {
  return x;
 }

 public double[][] getV() {
  return v;
 }

 public double[][] getF() {
  return f;
 }

 public double[][] getXwrap() {
  return xwrap;
 }

 public double getUs() {
  return Us;
 }

 public double getVir() {
  return Vir;
 }

 public double getP() {
  return p;
 }

 public double getUshift() {
  return ushift;
 }

 public double getUtail() {
  return Utail;
 }

 public double getPtail() {
  return ptail;
 }

 public int getStep() {
  return step;
 }

 public Ave getAvp() {
  return avp;
 }

 public double[] getXij() {
  return xij;
 }

 public Random getRng() {
  return rng;
 }

 public double getDen() {
  return den;
 }

 public void setDen(double den) {
  this.den = den;
 }
}

/** Average and standard deviation */
class Ave {
 double cnt, xsm, x2sm;

 public void clear() {
  cnt = xsm = x2sm = 0;
 }

 public void add(double x) {
  cnt += 1;
  xsm += x;
  x2sm += x * x;
 }

 public double getAve() {
  return (cnt > 0) ? xsm / cnt : 0;
 }

 public double getVar() {
  if (cnt <= 1)
   return 0;
  double av = xsm / cnt;
  return x2sm / cnt - av * av;
 }

 public double getStd() {
  return Math.sqrt(getVar());
 }
}

/** Panel for drawing particles */
class XYZCanvas extends JPanel implements MouseListener, MouseMotionListener, MouseWheelListener {
 Image img;
 Graphics imgG;
 Dimension imgSize;

 double realSpan; // size of a real span (saved as a copy)
 double zoomScale = 1.0;
 public Matrix3D viewMatrix = new Matrix3D(); // view matrix
 private Matrix3D tmpMatrix = new Matrix3D(); // temporary matrix
 int mouseX, mouseY; // mouse position
 XYZModel model;

 public static final long serialVersionUID = 2L;

 public XYZCanvas() {
  super();
  addMouseListener(this);
  addMouseMotionListener(this);
  addMouseWheelListener(this);
 }

 /**
  * Prepare a buffer for the image, return if the canvas is ready Only need
  * to call this when the size of the canvas is changed, Since we
  * automatically detect the size change in paint() only call this on start
  */
 public boolean newImgBuf() {
  Dimension sz = getSize();
  if (sz.width == 0 || sz.height == 0)
   return false;
  // quit if the current image already has the right size
  if (img != null && imgG != null && sz.equals(imgSize))
   return true;
  img = createImage(sz.width, sz.height);
  if (imgG != null)
   imgG.dispose();
  imgG = img.getGraphics();
  imgSize = sz;
  return true;
 }

 public void update(Graphics g) {
  if (img == null)
   g.clearRect(0, 0, getSize().width, getSize().height);
  paintComponent(g);
 }

 protected void paintComponent(Graphics g) {
  super.paintComponent(g);
  if (model != null) {
   newImgBuf(); // refresh the image buffer if necessary
   // compute the real-to-screen ratio, this variable differs
   // from model.real2Screen by zoomScale
   Dimension sz = getSize();
   double real2Screen0 = model.getScaleFromSpan(realSpan, sz.width, sz.height);
   model.setMatrix(viewMatrix, real2Screen0 * zoomScale, sz.width / 2, sz.height / 2);
   imgG.setColor(Color.BLACK);
   imgG.fillRect(0, 0, sz.width, sz.height);
   model.paint(imgG);
   g.drawImage(img, 0, 0, this);
  }
 }

 /**
  * Refresh the coordinates x[][] is the wrapped coordinates `n' may be less
  * than x.length
  */
 public void refresh(double x[][], int n, boolean center, boolean adjScale) {
  if (model == null) {
   model = new XYZModel();
   adjScale = true;
  }
  model.updateXYZ(x, n, center);
  if (adjScale)
   realSpan = model.getSpan(x, n);
  repaint();
 }

 /** Event handling */
 public void mouseClicked(MouseEvent e) {
 }

 public void mousePressed(MouseEvent e) {
  mouseX = e.getX();
  mouseY = e.getY();
  // consume this event so that it will not be processed
  // in the default manner
  e.consume();
 }

 public void mouseReleased(MouseEvent e) {
 }

 public void mouseEntered(MouseEvent e) {
 }

 public void mouseExited(MouseEvent e) {
 }

 public void mouseDragged(MouseEvent e) {
  int x = e.getX();
  int y = e.getY();
  tmpMatrix.unit();
  tmpMatrix.xrot(360.0 * (mouseY - y) / getSize().height);
  tmpMatrix.yrot(360.0 * (x - mouseX) / getSize().width);
  viewMatrix.mult(tmpMatrix);
  repaint();
  mouseX = x;
  mouseY = y;
  e.consume();
 }

 public void mouseMoved(MouseEvent e) {
 }

 public void mouseWheelMoved(MouseWheelEvent e) {
  int notches = e.getWheelRotation();
  if ((zoomScale -= 0.05f * notches) < 0.09999f)
   zoomScale = 0.1f;
  repaint();
 }
}

/** A set of atoms with 3D coordinates */
class XYZModel {
 double realXYZ[][]; // 3D real coordinates [np][3]
 int screenXYZ[][]; // 3D screen coordinates in pixels [np][3]
      // only the first two dimensions are used in drawing
 int zOrder[]; // z-order, larger is closer to the viewer
 int np = -1;

 boolean transformed;
 // rotation/scaling/translation matrix for the conversion from realXYZ[] to
 // screenXYZ[]
 Matrix3D mat = new Matrix3D();
 double real2Screen = 50.0; // real size --> screen size
 double ballSize = 0.5; // ball size (radius) in terms of the real
       // coordinates
       // 0.5 for hard spheres

 Atom atomDefault = new Atom(0.1, 0.7, 0.1, 1.0, 0.5, 0.5, 1.0);
 Atom atoms[]; // for colors of atoms

 XYZModel() {
 }

 /** Set the color of a particular atom */
 void setAtom(int i, Atom atom) {
  if (i >= 0 && i < atoms.length)
   atoms[i] = atom;
 }

 /**
  * Refresh coordinates x[0..n-1][3] is a three-dimensional vector n can be
  * less than x.length
  */
 void updateXYZ(double x[][], int n, boolean center) {
  if (n != np) {
   // System.out.printf("XYZModel.updateXYZ n %d --> %d, (%g, %g, %g)
   // (%g, %g, %g)\n", np, n, x[0][0], x[0][1], x[0][2], x[n-1][0],
   // x[n-1][1], x[n-1][2]);
   np = n;
   realXYZ = new double[np][3];
   screenXYZ = new int[np][3];
   zOrder = new int[np];
   atoms = new Atom[np];
   for (int i = 0; i < np; i++)
    atoms[i] = atomDefault;
  }

  for (int d = 0; d < 3; d++) {
   double xc = 0;
   if (center) {
    for (int i = 0; i < np; i++)
     xc += x[i][d];
    xc /= np;
   }
   for (int i = 0; i < np; i++)
    realXYZ[i][d] = x[i][d] - xc;
  }

  transformed = false;
 }

 /**
  * Set the view matrix `s' is the scaling factor of translating real
  * coordinates to the screen coordinates (x0, y0) the screen coordinates of
  * the center
  */
 void setMatrix(Matrix3D viewMat, double s, double x0, double y0) {
  mat.unit();
  mat.mult(viewMat);
  mat.scale(s, s, s);
  real2Screen = s;
  mat.translate(x0, y0, 0);
  transformed = false;
 }

 /**
  * Get the span of the model `n' may be less than x.length
  */
 double getSpan(double x[][], int n) {
  int dim = x[0].length;
  double realSpan = 0, del, fw, fh;

  for (int d = 0; d < dim; d++) {
   double xmin = 1e30, xmax = -1e30;
   for (int i = 0; i < n; i++)
    if (x[i][d] < xmin)
     xmin = x[i][d];
    else if (x[i][d] > xmax)
     xmax = x[i][d];
   if ((del = xmax - xmin) > realSpan)
    realSpan = del;
  }
  return realSpan;
 }

 /**
  * Translate a given span, return the real-to-screen ratio `w' and `h' are
  * the width and height of the screen in pixels
  */
 double getScaleFromSpan(double realSpan, int w, int h) {
  realSpan += ballSize * 2; // add two radii
  double fw = w / realSpan;
  double fh = h / realSpan;
  double facShrink = 0.9; // shrink a bit for the margin
  return (fw < fh ? fw : fh) * facShrink;
 }

 /** Compute the Z-order */
 void getZOrder() {
  // transform the coordinates
  if (!transformed) {
   mat.transform(realXYZ, screenXYZ, np);
   transformed = true;
  }

  // bubble sort z-order
  // zOrder[0] is the fartherest from the viewer
  // zOrder[np - 1] is the nearest to the viewer
  for (int i = 0; i < np; i++)
   zOrder[i] = i;
  for (int i = 0; i < np; i++) {
   // find the particle with the smallest z
   int jm = i, k;
   for (int j = i + 1; j < np; j++)
    if (screenXYZ[zOrder[j]][2] < screenXYZ[zOrder[jm]][2])
     jm = j;
   if (jm != i) {
    k = zOrder[i];
    zOrder[i] = zOrder[jm];
    zOrder[jm] = k;
   }
  }
 }

 /** Draw atom */
 void drawAtom(Graphics g, int iz) {
  int i = zOrder[iz];
  Atom atom = atoms[i];
  if (atom == null)
   return;
  double zMin = screenXYZ[zOrder[0]][2]; // fartherest from the viewer
  double zMax = screenXYZ[zOrder[np - 1]][2]; // closest to the viewer
  int greyScale = Atom.nBalls - 1;
  if (zMin != zMax) // zMin == zMax means the two-dimensional case
   greyScale = (int) (Atom.nBalls * (screenXYZ[i][2] - zMin) / (zMax - zMin) - 1e-6);
  // the atom closest to the viewer has a greyScale of Atom.nBalls - 1
  // the atom fartherest from the viewer has a greyScale of 0
  double radius = ballSize * atom.relRadius * real2Screen;
  atom.paint(g, screenXYZ[i][0], screenXYZ[i][1], greyScale, radius);
 }

 /** Paint this model to the graphics context `g' */
 void paint(Graphics g) {
  if (realXYZ == null || np <= 0)
   return;
  getZOrder();
  for (int iz = 0; iz < np; iz++)
   drawAtom(g, iz);
 }
}

/** Atom: a ball */
class Atom {
 private final static int R = 120;
 private final static int hx = 45; // (hx, hy) is the offset of the spot
          // light from the center
 private final static int hy = 45;
 private static int maxr; // maximal distance from the spotlight
 public final static int nBalls = 16; // shades of grey
 // spotlight brightness (0.0, 1.0)
 // 1.0 means the spotlight is pure white
 // 0.0 means the spotlight is the same as the solid color
 double spotlightAmp = .4;
 // color damping along the distance from the spotlight
 double rContrast = 0.7;
 // z-depth contrast (0.0, inf)
 // inf means the fartherest atom is completely dark
 // 0.0 means the fartherest atom is the same as the
 // nearest atom
 double zContrast = 2.0;

 private static byte data[];

 static { // data[] is a bitmap image of the ball of radius R
  data = new byte[R * 2 * R * 2];
  for (int Y = -R; Y < R; Y++) {
   int x0 = (int) (Math.sqrt(R * R - Y * Y) + 0.5);
   for (int X = -x0; X < x0; X++) {
    // sqrt(x^2 + y^2) gives distance from the spot light
    int x = X + hx, y = Y + hy;
    int r = (int) (Math.sqrt(x * x + y * y) + 0.5);
    // set the maximal intensity to the maximal distance
    // (in pixels) from the spot light
    if (r > maxr)
     maxr = r;
    data[(Y + R) * (R * 2) + (X + R)] = (r <= 0) ? 1 : (byte) r;
   }
  }
 }

 // the following variables are atom dependent
 private int Rl = 100, Gl = 100, Bl = 100;
 private Image balls[]; // 0..nBalls-1, at different z distances

 /** Constructor */
 Atom(int r, int g, int b) {
  setRGB(r, g, b);
 }

 /** colors that range from 0 to 1 */
 Atom(double r, double g, double b) {
  setRGB(r, g, b);
 }

 double relRadius = 1; // only used to save the radius

 Atom(double r, double g, double b, double rad) {
  this(r, g, b);
  relRadius = rad;
 }

 Atom(double r, double g, double b, double rad, double rcontrast, double spotlight, double zcontrast) {
  relRadius = rad;
  rContrast = rcontrast;
  spotlightAmp = spotlight;
  zContrast = zcontrast;
  setRGB(r, g, b);
 }

 /** Set color */
 void setRGB(int r, int g, int b) {
  Rl = r;
  Gl = g;
  Bl = b;
  makeBalls();
 }

 void setRGB(double r, double g, double b) {
  Rl = (int) (256 * r - 1e-6);
  Gl = (int) (256 * g - 1e-6);
  Bl = (int) (256 * b - 1e-6);
  makeBalls();
 }

 /** Linearly interpolate colors */
 private int blend(double fg, double bg, double fgfactor) {
  return (int) (bg + (fg - bg) * fgfactor);
 }

 // need a component instance to call createImage()
 private static Component component = new Applet();

 /** Prepare ball images with different sizes */
 private void makeBalls() {
  balls = new Image[nBalls];
  byte red[] = new byte[256];
  byte green[] = new byte[256];
  byte blue[] = new byte[256];
  for (int id = 0; id < nBalls; id++) {
   // smaller `b' means closer to black
   // if id == 0 (fartherest from the viewer)
   // b = 1/(1 + zContrast)
   // the outer blend() gives a color close to bgGrey
   // if id == nBalls - 1 (closest to the viewer),
   // b = 1, the outer blend() gives the color of
   // the inner blend()
   double b = (zContrast * id / (nBalls - 1) + 1) / (zContrast + 1);
   for (int i = maxr; i >= 0; --i) {
    // closeness to the spotlight
    double q = 1 - 1. * i / maxr;
    // dampness of the color along the radius
    double p = 1 - rContrast * i / maxr;
    // contrast of the spotlight
    // if i == 0 (closest to the spotlight),
    // d = 1.0 - spotlightAmp, the inner
    // blend() gives a color close to 255
    // (if spotlightAmp == 1).
    // if i == maxr (fartherest from the spotlight),
    // d = 1.0, the inner blend() gives
    // the foreground color, i.e., Rl, Gl, Bl
    // Thus, the inner blend() depends on the distance
    // from the spotlight, i == 0 means to be closest
    // to the spotlight
    double d = 1 - q * spotlightAmp;
    red[i] = (byte) blend(blend(Rl * p, 255, d), 0, b);
    green[i] = (byte) blend(blend(Gl * p, 255, d), 0, b);
    blue[i] = (byte) blend(blend(Bl * p, 255, d), 0, b);
   }
   // 256 color model
   IndexColorModel model = new IndexColorModel(8, maxr + 1, red, green, blue, 0);
   balls[id] = component.createImage(new MemoryImageSource(R * 2, R * 2, model, data, 0, R * 2));
  }
 }

 /**
  * Draw a ball at screen coordinate (x, y) with a ball index `id' (0, 0)
  * represents the top-left corner `x', `y', `radius' are given in pixels the
  * ball index (gray code) `id' can be 0 to 15
  */
 void paint(Graphics gc, int x, int y, int id, double radius) {
  if (balls == null)
   makeBalls();
  Image img = balls[id]; // id = [0..15]

  int size = (int) (radius * 2 + .5);
  gc.drawImage(img, x - size / 2, y - size / 2, size, size, null);
  // System.out.println("" + x + " " + y + " " + id + " " + radius);
 }
}

class Matrix3D {
 double xx, xy, xz, xo;
 double yx, yy, yz, yo;
 double zx, zy, zz, zo;
 static final double pi = 3.14159265;

 /** Create a new unit matrix */
 Matrix3D() {
  xx = 1.0;
  yy = 1.0;
  zz = 1.0;
 }

 /** Scale along each axis independently */
 void scale(double xf, double yf, double zf) {
  xx *= xf;
  xy *= xf;
  xz *= xf;
  xo *= xf;
  yx *= yf;
  yy *= yf;
  yz *= yf;
  yo *= yf;
  zx *= zf;
  zy *= zf;
  zz *= zf;
  zo *= zf;
 }

 /** Translate the origin */
 void translate(double x, double y, double z) {
  xo += x;
  yo += y;
  zo += z;
 }

 /** Rotate theta degrees around the y axis */
 void yrot(double theta) {
  theta *= (pi / 180);
  double ct = Math.cos(theta);
  double st = Math.sin(theta);

  double Nxx = (xx * ct + zx * st);
  double Nxy = (xy * ct + zy * st);
  double Nxz = (xz * ct + zz * st);
  double Nxo = (xo * ct + zo * st);

  double Nzx = (zx * ct - xx * st);
  double Nzy = (zy * ct - xy * st);
  double Nzz = (zz * ct - xz * st);
  double Nzo = (zo * ct - xo * st);

  xo = Nxo;
  xx = Nxx;
  xy = Nxy;
  xz = Nxz;
  zo = Nzo;
  zx = Nzx;
  zy = Nzy;
  zz = Nzz;
 }

 /** Rotate theta degrees about the x axis */
 void xrot(double theta) {
  theta *= (pi / 180);
  double ct = Math.cos(theta);
  double st = Math.sin(theta);

  double Nyx = (yx * ct + zx * st);
  double Nyy = (yy * ct + zy * st);
  double Nyz = (yz * ct + zz * st);
  double Nyo = (yo * ct + zo * st);

  double Nzx = (zx * ct - yx * st);
  double Nzy = (zy * ct - yy * st);
  double Nzz = (zz * ct - yz * st);
  double Nzo = (zo * ct - yo * st);

  yo = Nyo;
  yx = Nyx;
  yy = Nyy;
  yz = Nyz;
  zo = Nzo;
  zx = Nzx;
  zy = Nzy;
  zz = Nzz;
 }

 /** Rotate theta degrees about the z axis */
 void zrot(double theta) {
  theta *= pi / 180;
  double ct = Math.cos(theta);
  double st = Math.sin(theta);

  double Nyx = (yx * ct + xx * st);
  double Nyy = (yy * ct + xy * st);
  double Nyz = (yz * ct + xz * st);
  double Nyo = (yo * ct + xo * st);

  double Nxx = (xx * ct - yx * st);
  double Nxy = (xy * ct - yy * st);
  double Nxz = (xz * ct - yz * st);
  double Nxo = (xo * ct - yo * st);

  yo = Nyo;
  yx = Nyx;
  yy = Nyy;
  yz = Nyz;
  xo = Nxo;
  xx = Nxx;
  xy = Nxy;
  xz = Nxz;
 }

 /** Multiply this matrix by a second: M = M*R */
 void mult(Matrix3D rhs) {
  double lxx = xx * rhs.xx + yx * rhs.xy + zx * rhs.xz;
  double lxy = xy * rhs.xx + yy * rhs.xy + zy * rhs.xz;
  double lxz = xz * rhs.xx + yz * rhs.xy + zz * rhs.xz;
  double lxo = xo * rhs.xx + yo * rhs.xy + zo * rhs.xz + rhs.xo;

  double lyx = xx * rhs.yx + yx * rhs.yy + zx * rhs.yz;
  double lyy = xy * rhs.yx + yy * rhs.yy + zy * rhs.yz;
  double lyz = xz * rhs.yx + yz * rhs.yy + zz * rhs.yz;
  double lyo = xo * rhs.yx + yo * rhs.yy + zo * rhs.yz + rhs.yo;

  double lzx = xx * rhs.zx + yx * rhs.zy + zx * rhs.zz;
  double lzy = xy * rhs.zx + yy * rhs.zy + zy * rhs.zz;
  double lzz = xz * rhs.zx + yz * rhs.zy + zz * rhs.zz;
  double lzo = xo * rhs.zx + yo * rhs.zy + zo * rhs.zz + rhs.zo;

  xx = lxx;
  xy = lxy;
  xz = lxz;
  xo = lxo;
  yx = lyx;
  yy = lyy;
  yz = lyz;
  yo = lyo;
  zx = lzx;
  zy = lzy;
  zz = lzz;
  zo = lzo;
 }

 /** Recover the unit matrix */
 void unit() {
  xo = 0;
  xx = 1;
  xy = 0;
  xz = 0;
  yo = 0;
  yx = 0;
  yy = 1;
  yz = 0;
  zo = 0;
  zx = 0;
  zy = 0;
  zz = 1;
 }

 /**
  * Transform np points from v into tv. v contains the input coordinates in
  * floating point. Three successive entries in the array constitute a point.
  * tv ends up holding the transformed points as integers; three successive
  * entries per point
  */
 void transform(double v[][], int tv[][], int np) {
  // np can be different from v.length
  for (int i = 0; i < np; i++) {
   double x = v[i][0], y = v[i][1], z = v[i][2];
   tv[i][0] = (int) (xx * x + xy * y + xz * z + xo);
   tv[i][1] = (int) (yx * x + yy * y + yz * z + yo);
   tv[i][2] = (int) (zx * x + zy * y + zz * z + zo);
  }
 }
}

2) EvolutionLibrary.java:

package md_3d;

/* Developed by Samuel Graván 2015 */

import java.awt.Graphics;
import java.util.ArrayList;
import java.util.Collections;
import java.util.Comparator;

import javax.swing.JPanel;

class Util {
 public static float Round(float Rval, int Rpl) {
  float p = (float) Math.pow(10, Rpl);
  Rval = Rval * p;
  float tmp = Math.round(Rval);
  return (float) tmp / p;
 }

}

class ValorComparator implements Comparator {

 public int compare(Object o1, Object o2) {
  Individuo a = (Individuo) o1;
  Individuo b = (Individuo) o2;

  //return Double.compare(a.getFreedom(), b.getFreedom());
  return Double.compare(a.getFreedom(), b.getFreedom());
 }

 public boolean equals(Object o) {
  return this == o;
 }
}

class Individuo {
 private LJ3MD lj;
 private double value;
 private double freedom;

 public double getFreedom() {
  return freedom;
 }

 public void setFreedom(double freedom) {
  this.freedom = freedom;
 }

 public void setValue(double value) {
  this.value = value;
 }

 Individuo(LJ3MD lj) {
  this.lj = lj;
  value = 0d;
  freedom = 0d;
 }

 public double getValue() {
  return value;
 }

 public LJ3MD getLj() {
  return lj;
 }

 public void setLj(LJ3MD lj) {
  this.lj = lj;
 }
}

class Poblacion {
 private java.util.List<Individuo> individuos;

 Poblacion() {
  individuos = new ArrayList<Individuo>();
 }

 Poblacion(long i, LJ3MD mdI) {
  individuos = new ArrayList<Individuo>();

  for (long j = 0; j < i; j++) {
   LJ3MD md = new LJ3MD(0);

   md.init(mdI.getDen(), mdI.getN(), mdI.getV(), mdI.getF(), mdI.getX(), mdI.getStep(), mdI.getRho(),
     mdI.getDt(), mdI.getVol(), mdI.getRc(), mdI.getkT(), mdI.getXWrap(), mdI.getK(), mdI.getU(),
     mdI.getUs(), mdI.getE(), mdI.getVir(), mdI.getP(), mdI.getUshift(), mdI.getPtail(), mdI.getUtail(),
     mdI.getAvU(), mdI.getAvK(), mdI.getAvp(), mdI.getXij());

   Individuo I = new Individuo(md);
   I.setFreedom(mdI.getMeanFreedom());
   I.setValue(mdI.getMeanEnergy());

   individuos.add(I);
  }
 }

 public void evaluar(int steps) {

  for (Individuo i : this.getIndividuos()) {

   double Ei = i.getLj().getE();

   for (int k = 0; k <= steps; k++) {
    i.getLj().vv();
   }

   double Ej = i.getLj().getE();

   double BQ = (-1) * (Ei - Ej);

   i.setValue(i.getValue() + BQ);
   i.setFreedom((i.getFreedom() + i.getLj().getMeanFreedom()) / steps);
  }
 }

 public void nuevoIndividuo(Poblacion p) {
  this.individuos.addAll(p.getIndividuos());
 }

 public void nuevoIndividuo(java.util.List<Individuo> p) {
  this.individuos.addAll(p);
 }

 public java.util.List<Individuo> getIndividuos() {
  return this.individuos;
 }

 public void setPoblacion(java.util.List<Individuo> individuos) {
  this.individuos = individuos;
 }
}

public class EvolutionLibrary {

 private long N = 250;
 private long g = 1250; // 500 generaciones
 private int steps = 100; // 1 segundo
 Poblacion P;
 
 public EvolutionLibrary(long N, long g, int steps) {
  this.N = N;
  this.g = g;
  this.steps = steps;
 }

 public java.util.List<Individuo> mutacion() {
  java.util.List<Individuo> mutaciones = new ArrayList<Individuo>();
  for (Individuo i : P.getIndividuos()) {
   LJ3MD md = new LJ3MD(0);

   md.init(i.getLj().getDen(), i.getLj().getN(), i.getLj().getV(), i.getLj().getF(), i.getLj().getX(),
     i.getLj().getStep(), i.getLj().getRho(), i.getLj().getDt(), i.getLj().getVol(), i.getLj().getRc(),
     i.getLj().getkT(), i.getLj().getXWrap(), i.getLj().getK(), i.getLj().getU(), i.getLj().getUs(),
     i.getLj().getE(), i.getLj().getVir(), i.getLj().getP(), i.getLj().getUshift(), i.getLj().getPtail(),
     i.getLj().getUtail(), i.getLj().getAvU(), i.getLj().getAvK(), i.getLj().getAvp(),
     i.getLj().getXij());

   md.vvEvolution();

   Individuo I = new Individuo(md);
   I.setFreedom(i.getFreedom());
   I.setValue(i.getValue());

   mutaciones.add(I);
  }

  return mutaciones;
 }

 public Poblacion buscarSelection(LJ3MD mdI, XYZCanvas canvas, LJ3MDApp app, Graphics gG) {
  int t = 0;

  P = new Poblacion(N, mdI);

  P.evaluar(steps);
  do {

   Poblacion Q = new Poblacion();
   Q.setPoblacion(mutacion());

   P.nuevoIndividuo(Q);

   P.evaluar(steps);

   java.util.List<Individuo> resultados = new ArrayList<Individuo>();
   java.util.List<Individuo> individuos = new ArrayList<Individuo>(P.getIndividuos());

   ValorComparator vc = new ValorComparator();
   Collections.sort(individuos, vc);
   long j = 0;

   for (Individuo i : individuos) {
    if (j < individuos.size() / 2) {
     resultados.add(i);
     j++;
    } else
     break;
   }

   P.setPoblacion(resultados);
   
   // Refrescamos la pantalla principal con el mejor individuo de esta generación  
   app.md = new LJ3MD(0);
   app.md.init(resultados.get(0).getLj().getDen(), resultados.get(0).getLj().getN(), resultados.get(0).getLj().getV(), resultados.get(0).getLj().getF(), resultados.get(0).getLj().getX(),
     resultados.get(0).getLj().getStep(), resultados.get(0).getLj().getRho(), resultados.get(0).getLj().getDt(), resultados.get(0).getLj().getVol(), resultados.get(0).getLj().getRc(),
     resultados.get(0).getLj().getkT(), resultados.get(0).getLj().getXWrap(), resultados.get(0).getLj().getK(), resultados.get(0).getLj().getU(), resultados.get(0).getLj().getUs(),
     resultados.get(0).getLj().getE(), resultados.get(0).getLj().getVir(), resultados.get(0).getLj().getP(), resultados.get(0).getLj().getUshift(), resultados.get(0).getLj().getPtail(),
     resultados.get(0).getLj().getUtail(), resultados.get(0).getLj().getAvU(), resultados.get(0).getLj().getAvK(), resultados.get(0).getLj().getAvp(),
     resultados.get(0).getLj().getXij());
   
   canvas.refresh(resultados.get(0).getLj().getXWrap(), resultados.get(0).getLj().N, true, false);
   canvas.paint(gG);
   app.update(gG);
   app.repaint();
   app.getContentPane().update(gG);
   
   t++;
  } while (t < g);

  return P;
 }
}