```python
import numpy as np
import matplotlib.pyplot as plt
import sympy as sy
import simtk.unit as unit
from simtk import openmm as mm
from simtk.openmm import app
import skopt as skopt
from tqdm import tqdm
```

# A Lennard-Jones Fluid

## The Lennard-Jones potential

The Lennard-Jones (LJ) potential between two particles is defined by the following equation, where $x$ is the distance between the particles, and $\sigma$ and $\epsilon$ are two parameters of the potential:
    
\begin{equation}
V(x) = 4 \epsilon \left[ \left( \frac{\sigma}{x} \right)^{12} - \left( \frac{\sigma}{x} \right)^6 \right]
\end{equation}

Lets see the shape of this function:


```python
def LJ (x, sigma, epsilon):
    
    t = sigma/x
    t6 = t**6
    t12 = t6**2
    
    return 4.0*epsilon*(t12-t6)
```


```python
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

xlim_figure = [0.01, 6.0]
ylim_figure = [-2.0, 10.0]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstrom
plt.plot(x, LJ(x, sigma, epsilon))
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
```

The way the LJ potential is built, the $\sigma$ and $\epsilon$ parameters have a straightforward interpretation. The cut with $y=0$ is located in $x=\sigma$:


```python
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

xlim_figure = [0.01, 6.0]
ylim_figure = [-2.0, 10.0]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstrom
plt.plot(x, LJ(x, sigma, epsilon))
plt.hlines(0, xlim_figure[0], xlim_figure[1], linestyles='dotted', color='gray')
plt.vlines(sigma._value, ylim_figure[0], ylim_figure[1], linestyles='dashed', color='red')
plt.text(sigma._value+0.02*xlim_figure[1], 0.7*ylim_figure[1], '$\sigma$', fontsize=14)
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
```

And $\epsilon$ is the depth of the minimum measured from $y=0$:


```python
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

xlim_figure = [0.01, 6.0]
ylim_figure = [-2.0, 10.0]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstrom
plt.plot(x, LJ(x, sigma, epsilon))
plt.hlines(0, xlim_figure[0], xlim_figure[1], linestyles='dotted', color='gray')
plt.hlines(-epsilon._value, xlim_figure[0], xlim_figure[1], linestyles='dashed', color='red')
plt.annotate(text='', xy=(1.0,0.0), xytext=(1.0,-epsilon._value), arrowprops=dict(arrowstyle='<->'))
plt.text(1.0+0.02*xlim_figure[1], -0.7*epsilon._value, '$\epsilon$', fontsize=14)
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
```

Notice that the LJ potential has physical meaning when $\epsilon>0$ and $\sigma>0$ only. Actually, the potential vanishes whether $\epsilon=0$ or $\sigma=0$.

### The Lennard Jones minimum and the size of the particles

The LJ potential has a single minimum located in $x_{min}$. Lets equal to $0$ the first derivative of the potential to find the value of $x_{min}$:


```python
x, sigma, epsilon = sy.symbols('x sigma epsilon', real=True, positive=True)
V = 4.0*epsilon*((sigma/x)**12-(sigma/x)**6)
gradV = sy.diff(V,x)
roots=sy.solve(gradV, x)
x_min = roots[0]
```


```python
x_min
```




$\displaystyle 1.12246204830937 \sigma$



The minimum is then located in:

\begin{equation}
x_{min} = 2^{1/6} \sigma
\end{equation}

where the potential takes the value:

\begin{equation}
V(x_{min}) = -\epsilon
\end{equation}


```python
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

x_min = 2**(1/6)*sigma
y_min = -epsilon

xlim_figure = [x_min._value-0.4, x_min._value+0.4]
ylim_figure = [y_min._value-0.1, y_min._value+0.5]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstroms
plt.plot(x, LJ(x, sigma, epsilon))
plt.hlines(y_min._value, xlim_figure[0], xlim_figure[1], linestyles='dashed', color='gray')
plt.vlines(x_min._value, ylim_figure[0], ylim_figure[1], linestyles='dashed', color='gray')
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
```

This way two particles in the equilibrium position will be placed at a $2^{1/6} \sigma$ distance. The potential is thereby modeling two "soft spheres" atracting each other very lightly. Their radii, given that both particles are equal, are equal to $r$:

\begin{equation}
r = \frac{1}{2} x_{min} = 2^{-5/6} \sigma
\end{equation}

And we say these spheres are "soft" because their volume is not limited by a hard-wall potential, they can penetrate each other suffering a not infinite repulsive force.

### Time period of the small harmonic oscillations around the minimum

If we want to perform a molecular simulation of this two particles we should wonder how big the integrator timestep must be. To answer this question we can study the harmonic approximation around the minimum. Lets calculate the time period, $\tau$, of a small harmonic oscillation around the minimum:


```python
x, sigma, epsilon = sy.symbols('x sigma epsilon', real=True, positive=True)
V = 4.0*epsilon*((sigma/x)**12-(sigma/x)**6)
gradV = sy.diff(V,x)
grad2V = sy.diff(V,x,x)

x_min = sy.solve(gradV,x)[0]
k_harm = grad2V.subs(x, x_min)
```


```python
k_harm
```




$\displaystyle \frac{57.1464378708551 \epsilon}{\sigma^{2}}$



The harmonic constant of the second degree Taylor polynomial of the LJ potential at $x=x_{min}$ is then:

\begin{equation}
k_{harm} = 36·2^{2/3} \frac{\epsilon}{\sigma^2}
\end{equation}

The oscillation period of a particle with $m$ mass in an harmonic potential defined by $\frac{1}{2} k x²$ is:

\begin{equation}
\tau = 2 \pi \sqrt{ \frac{m}{k}}
\end{equation}

As such, the period of the small harmonic oscillations around the LJ minimum of particle with $m$ mass is:

\begin{equation}
\tau = 2 \pi \sqrt{ \frac{m}{k_{harm}}} = \frac{\pi}{3·2^{1/3}} \sqrt{\frac{m\sigma^2}{\epsilon}}
\end{equation}

With the mass and parameters taking values of amus, angstroms and kilocalories per mole, the time period is in the order of:


```python
mass = 50.0 * unit.amu
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

k = 36 * 2**(2/3) * epsilon/sigma**2

tau = 2*np.pi * np.sqrt(mass/k)

print(tau)
```

    0.5746513694274475 ps


But, is this characteristic time a good threshold for a LJ potential? If the oscillations around the minimum are not small enough, the harmonic potential of the second degree term of the taylor expansion is easily overcome by the sharp left branch of the LJ potential:


```python
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

k = 36 * 2**(2/3) * epsilon/sigma**2

x_min = 2**(1/6)*sigma
y_min = -epsilon

xlim_figure = [x_min._value-0.2, x_min._value+0.2]
ylim_figure = [y_min._value-0.1, y_min._value+0.6]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstroms
plt.plot(x, LJ(x, sigma, epsilon))
plt.plot(x, 0.5*k*(x-x_min)**2+y_min)
plt.hlines(y_min._value, xlim_figure[0], xlim_figure[1], linestyles='dashed', color='gray')
plt.vlines(x_min._value, ylim_figure[0], ylim_figure[1], linestyles='dashed', color='gray')
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
```

Let's imagine the following situation. Let a particle be in the harmonic potential at temperature of 300K. Will the particle be more constrained in space than in the well of the LJ potential? Will the particle feel the harmonic potential softer or sharper than the LJ? Lets make some numbers to evaluate if the oscillation time period of the harmonic approximation can be a good time threshold for the integration timestep of a molecular dynamics of the LJ potential.

The standard deviation of an harmonic oscillation with the shape $\frac{1}{2}k x^2$ in contact with a stochastic thermal bath can be computed as:

\begin{equation}
\beta = \frac{1}{k_{\rm B} T} 
\end{equation}

\begin{equation}
Z_x = \int_{-\infty}^{\infty} {\rm e}^{- \beta \frac{1}{2}k x^2} = \sqrt{\frac{2 \pi}{\beta k}}
\end{equation}

\begin{equation}
\left< x \right> = \frac{1}{Z_x} \int_{-\infty}^{\infty} x {\rm e}^{-\beta \frac{1}{2}k x^2} = 0
\end{equation}

\begin{equation}
\left< x^2 \right> = \frac{1}{Z_x} \int_{-\infty}^{\infty} x^{2} {\rm e}^{-\beta \frac{1}{2}k x^2} = \frac{1}{Z_x}  \sqrt{\frac{2 \pi}{\beta³ k^3}} = \frac{1}{\beta k}
\end{equation}


\begin{equation}
{\rm std} = \left( \left< x^2 \right> -\left< x \right>^2 \right)^{1/2} = \sqrt{ \frac{k_{\rm B}T}{k} }
\end{equation}


This way, in the case of the harmonic potential obtained as the second degree term of the Taylor expansion around the LJ minimum:


```python
mass = 50.0 * unit.amu
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole
temperature = 300 * unit.kelvin
kB = unit.BOLTZMANN_CONSTANT_kB * unit.AVOGADRO_CONSTANT_NA

k = 36 * 2**(2/3) * epsilon/sigma**2
std = np.sqrt(kB*temperature/k)

x_min = 2**(1/6)*sigma
y_min = -epsilon

xlim_figure = [x_min._value-0.4, x_min._value+0.4]
ylim_figure = [y_min._value-0.1, y_min._value+0.6]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstroms
plt.plot(x, LJ(x, sigma, epsilon))
plt.plot(x, 0.5*k*(x-x_min)**2+y_min)
plt.hlines(y_min._value, xlim_figure[0], xlim_figure[1], linestyles='dashed', color='gray')
plt.vlines(x_min._value, ylim_figure[0], ylim_figure[1], linestyles='dashed', color='gray')
plt.axvspan(x_min._value - std._value, x_min._value + std._value, alpha=0.2, color='red')
plt.annotate(text='', xy=(x_min._value, y_min._value - 0.5*(y_min._value-ylim_figure[0])),
             xytext=(x_min._value-std._value, y_min._value - 0.5*(y_min._value-ylim_figure[0])),
             arrowprops=dict(arrowstyle='<->'))
plt.text(x_min._value-0.6*std._value, y_min._value - 0.4*(y_min._value-ylim_figure[0]), '$std$', fontsize=14)
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
```

The harmonic potential is too soft as approximation. Its oscillation time used as threshold to choose the integration timestep can yield to numeric problems. Let's try with a stiffer potential, let's double the harmonic constant:


```python
mass = 50.0 * unit.amu
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole
temperature = 300 * unit.kelvin
kB = unit.BOLTZMANN_CONSTANT_kB * unit.AVOGADRO_CONSTANT_NA

k = 36 * 2**(2/3) * epsilon/sigma**2
std = np.sqrt(kB*temperature/k)

x_min = 2**(1/6)*sigma
y_min = -epsilon

xlim_figure = [x_min._value-0.4, x_min._value+0.4]
ylim_figure = [y_min._value-0.1, y_min._value+0.6]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstroms
plt.plot(x, LJ(x, sigma, epsilon))
plt.plot(x, 0.5*k*(x-x_min)**2+y_min)
plt.plot(x, k*(x-x_min)**2+y_min, label='2k_{harm}')
plt.hlines(y_min._value, xlim_figure[0], xlim_figure[1], linestyles='dashed', color='gray')
plt.vlines(x_min._value, ylim_figure[0], ylim_figure[1], linestyles='dashed', color='gray')
plt.axvspan(x_min._value - std._value, x_min._value + std._value, alpha=0.2, color='red')
plt.annotate(text='', xy=(x_min._value, y_min._value - 0.5*(y_min._value-ylim_figure[0])),
             xytext=(x_min._value-std._value, y_min._value - 0.5*(y_min._value-ylim_figure[0])),
             arrowprops=dict(arrowstyle='<->'))
plt.text(x_min._value-0.6*std._value, y_min._value - 0.4*(y_min._value-ylim_figure[0]), '$std$', fontsize=14)
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
```

Lets take then, as reference, an harmonic potential with constant equal to $2k_{harm}$ could be a better idea. Lets compute then the new time threshold to choose the integration timestep:

\begin{equation}
\tau = 2 \pi \sqrt{ \frac{m}{2k_{harm}}} = \frac{\pi}{3·2^{5/6}} \sqrt{\frac{m\sigma^2}{\epsilon}}
\end{equation}


```python
mass = 50.0 * unit.amu
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

k = 36 * 2**(2/3) * epsilon/sigma**2

tau = 2*np.pi * np.sqrt(mass/(2*k))

print(tau)
```

    0.4063398801402841 ps


It is an accepted rule of thumb that the integration timestep must be as large as $\tau / 10$, being $\tau$ the oscillation time period of the fastest possible vibration mode. So finally, in this case the integration time step should not be longer than:


```python
mass = 50.0 * unit.amu
sigma = 2.0 * unit.angstrom
epsilon = 1.0 * unit.kilocalories_per_mole

k = 36 * 2**(2/3) * epsilon/sigma**2

tau = 2*np.pi * np.sqrt(mass/(2*k))

print(tau/10.0)
```

    0.04063398801402841 ps


For the case of a LJ potential modeling an Argon fluid:


```python
mass = 39.9 * 2 * unit.amu # reduced mass of a pair of atoms
sigma = 3.4 * unit.angstrom
epsilon = 0.238 * unit.kilocalories_per_mole

k = 36 * 2**(2/3) * epsilon/sigma**2

tau = 2*np.pi * np.sqrt(mass/(2*k))

print(tau/10.0)
```

    0.17888187339607967 ps


## Two Lennard-Jones atoms in vacuum

### Ar-Ar


```python
mass = 39.948 * unit.amu
sigma = 3.404 * unit.angstroms
epsilon = 0.238 * unit.kilocalories_per_mole
charge = 0.0 * unit.elementary_charge
```


```python
system = mm.System()

non_bonded_force = mm.NonbondedForce()
non_bonded_force.setNonbondedMethod(mm.NonbondedForce.NoCutoff)

# First Ar atom
system.addParticle(mass)
non_bonded_force.addParticle(charge, sigma, epsilon)

# Second Ar atom
system.addParticle(mass)
non_bonded_force.addParticle(charge, sigma, epsilon)

system.addForce(non_bonded_force)

integrator = mm.VerletIntegrator(2*unit.femtoseconds)
platform = mm.Platform.getPlatformByName('CUDA')

context = mm.Context(system, integrator, platform)
```


```python
positions = np.zeros([2,3], float) * unit.angstroms

x = np.linspace(1.0, 8.0, 200, endpoint=True) * unit.angstroms
V = [] * unit.kilocalories_per_mole

for xi in x:
    positions[1,0] = xi
    context.setPositions(positions)
    state = context.getState(getEnergy=True)
    potential_energy = state.getPotentialEnergy()
    V.append(potential_energy)
```

Position of minimum:

\begin{equation}
x_{min} = 2^{1/6} \sigma
\end{equation}


```python
x_min = 2**(1/6)*sigma
```


```python
x_min
```




    Quantity(value=3.8208608124451056, unit=angstrom)




```python
k = 36 * 2**(2/3) * epsilon/sigma**2
reduced_mass = (mass**2)/(2*mass)

tau = 2*np.pi * np.sqrt(reduced_mass/k)

print(tau)
```

    1.267135457848527 ps



```python
V._value = np.array(V._value)

xlim_figure = [0.01, 8.0]
ylim_figure = [-2.0, 10.0]

plt.plot(x, V, linewidth=2)
plt.plot(x, LJ(x, sigma, epsilon), linestyle='--', color='red')
plt.vlines(x_min._value, ylim_figure[0], ylim_figure[1], linestyles='dashed', color='gray')
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
```


```python
xlim_figure = [3.0, 5.0]
ylim_figure = [-0.4, 0.2]

x = np.linspace(xlim_figure[0], xlim_figure[1], 100, True) * unit.angstrom

plt.plot(x, LJ(x, sigma, epsilon))
plt.vlines(x_min._value, ylim_figure[0], ylim_figure[1], linestyles='dashed', color='gray')

plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
```


```python
initial_distance = x_min + 0.05 * unit.angstroms
simulation_time = 4*tau
saving_time = 0.01*tau
integration_timestep = 0.01*tau

saving_steps = int(saving_time/integration_timestep)
n_savings = int(simulation_time/saving_time)

integrator = context.getIntegrator()
integrator.setStepSize(integration_timestep)

trajectory = []*unit.angstroms
times = []*unit.picoseconds

positions = np.zeros([2,3], float) * unit.angstroms
positions[1,0] = initial_distance
context.setPositions(positions)

velocities = np.zeros([2,3], float) * unit.angstroms/unit.picoseconds
context.setVelocities(velocities)

context.setTime(0.0*unit.picoseconds)

state = context.getState(getPositions=True)
positions = state.getPositions(asNumpy=True)
time = state.getTime()
distance = positions[1,0] - positions[0,0]
trajectory.append(distance)
times.append(time)

for _ in range(n_savings):
    
    integrator.step(saving_steps)
    state = context.getState(getPositions=True)
    positions = state.getPositions(asNumpy=True)
    time = state.getTime()
    distance = positions[1,0] - positions[0,0]
    trajectory.append(distance)
    times.append(time)

times._value = np.array(times._value)
trajectory._value = np.array(trajectory._value)

plt.plot(times, trajectory)
plt.show()
```

- Try to place the initial distance between the pair of atoms further than the simulated above. 

### Ar-Xe molecular system


```python
mass_Ar = 39.948 * unit.amu
sigma_Ar = 3.404 * unit.angstroms
epsilon_Ar = 0.238 * unit.kilocalories_per_mole
charge_Ar = 0.0 * unit.elementary_charge

mass_Xe = 131.293 * unit.amu
sigma_Xe = 3.961 * unit.angstroms
epsilon_Xe = 0.459 * unit.kilocalories_per_mole
charge_Xe = 0.0 * unit.elementary_charge
```


```python
system = mm.System()

non_bonded_force = mm.NonbondedForce()
non_bonded_force.setNonbondedMethod(mm.NonbondedForce.NoCutoff)

# Ar atom
system.addParticle(mass_Ar)
non_bonded_force.addParticle(charge_Ar, sigma_Ar, epsilon_Ar)

# Xe atom
system.addParticle(mass_Xe)
non_bonded_force.addParticle(charge_Xe, sigma_Xe, epsilon_Xe)

system.addForce(non_bonded_force)

integrator = mm.VerletIntegrator(2*unit.femtoseconds)
platform = mm.Platform.getPlatformByName('CUDA')

context = mm.Context(system, integrator, platform)
```


```python
positions = np.zeros([2,3], float) * unit.angstroms

x = np.linspace(1.0, 8.0, 200, endpoint=True) * unit.angstroms
V = [] * unit.kilocalories_per_mole

for xi in x:

    positions[1,0] = xi
    context.setPositions(positions)
    state = context.getState(getEnergy=True)
    potential_energy = state.getPotentialEnergy()
    V.append(potential_energy)
```


```python
sigma = (sigma_Ar+sigma_Xe)/2.0
epsilon = np.sqrt(epsilon_Ar*epsilon_Xe)

V._value = np.array(V._value)

xlim_figure = [0.01, 8.0]
ylim_figure = [-2.0, 10.0]

plt.plot(x, V, linewidth=2)
plt.plot(x, LJ(x, sigma, epsilon), linestyle='--', color='red')
plt.xlim(xlim_figure)
plt.ylim(ylim_figure)
plt.xlabel('x [{}]'.format(x.unit.get_symbol()))
plt.ylabel('V [{}]'.format(epsilon.unit.get_symbol()))
plt.show()
```

## Liquid Argon model


```python
n_particles = 1000
reduced_density = 0.75
mass = 39.948 * unit.amu
sigma = 3.404 * unit.angstroms
epsilon = 0.238 * unit.kilocalories_per_mole
charge = 0.0 * unit.elementary_charge

temperature = 300.0 * unit.kelvin

integration_timestep = 2.0 * unit.femtoseconds
collisions_rate = 1.0 / unit.picoseconds

equilibration_time = 1.0 * unit.nanoseconds

production_time = 5.0 * unit.nanoseconds
saving_time = 50.0 * unit.picoseconds
```

The Van der Waals radius of Argon is 1.88 angstroms.


```python
radius = 2.0**(-5/6) * sigma
print(radius)
```


```python
volume_particles = n_particles * sigma**3
volume = volume_particles/reduced_density
l_box = volume**(1/3)
```


```python
system = mm.System()
```


```python
v1 = np.zeros(3) * unit.angstroms
v2 = np.zeros(3) * unit.angstroms
v3 = np.zeros(3) * unit.angstroms

v1[0] = l_box
v2[1] = l_box
v3[2] = l_box

system.setDefaultPeriodicBoxVectors(v1, v2, v3)
```


```python
non_bonded_force = mm.NonbondedForce()
non_bonded_force.setNonbondedMethod(mm.NonbondedForce.CutoffPeriodic)
non_bonded_force.setCutoffDistance(3.0*sigma)
non_bonded_force.setUseSwitchingFunction(True)
non_bonded_force.setSwitchingDistance(2.0*sigma)
non_bonded_force.setUseDispersionCorrection(True)
```


```python
for _ in range(n_particles):
    system.addParticle(mass)
    non_bonded_force.addParticle(charge, sigma, epsilon)
```


```python
_ = system.addForce(non_bonded_force)
```


```python
space = skopt.Space([[0.0, l_box._value], [0.0, l_box._value], [0.0, l_box._value]])
```


```python
integrator = mm.LangevinIntegrator(temperature, collisions_rate, integration_timestep)
platform = mm.Platform.getPlatformByName('CUDA')
context = mm.Context(system, integrator, platform)
```


```python
grid_generator = skopt.sampler.Grid(use_full_layout=False)
initial_positions = grid_generator.generate(space.dimensions, n_particles)
initial_positions = np.array(initial_positions)*unit.angstroms
```


```python
context.setPositions(initial_positions)
context.setVelocitiesToTemperature(temperature)
```


```python
state=context.getState(getEnergy=True)
print("Before minimization: {}".format(state.getPotentialEnergy()))
mm.LocalEnergyMinimizer_minimize(context)
state=context.getState(getEnergy=True)
print("After minimization: {}".format(state.getPotentialEnergy()))
```


```python
equilibration_n_steps = int(equilibration_time/integration_timestep)
integrator.step(equilibration_n_steps)
context.setTime(0.0*unit.picoseconds)
```


```python
production_n_steps = int(production_time/integration_timestep)
saving_n_steps = int(saving_time/integration_timestep)
n_saving_periods = int(production_n_steps/saving_n_steps)

time = np.zeros([n_saving_periods]) * unit.nanoseconds
trajectory = np.zeros([n_saving_periods, n_particles, 3]) * unit.angstroms
potential_energy = np.zeros([n_saving_periods]) * unit.kilocalories_per_mole

for ii in tqdm(range(n_saving_periods)):
    integrator.step(saving_n_steps)
    state = context.getState(getPositions=True, getEnergy=True)
    time[ii] = state.getTime()
    trajectory[ii,:,:] = state.getPositions(asNumpy=True)
    potential_energy = state.getPotentialEnergy()
```


```python
trajectory_mem = trajectory.size * trajectory.itemsize * unit.bytes
print('Trajectory size: {} GB'.format(trajectory_mem._value/(1024*1024)))
```


```python
l_box
```


```python
trajectory.max()
```


```python
- Calcular difusion de las partículas
- Calcular RDF
- Energy
```

0456b4d03eb0d0aa9ca356d326da11587502340f
ipynb
Jupyter Notebook
docs/contents/LJ_simulation.ipynb
uibcdf/OpenMembrane
c9705cb32706b882bdb3d75d19539ca323e6b741
MIT
docs/contents/LJ_simulation.ipynb
uibcdf/OpenMembrane
c9705cb32706b882bdb3d75d19539ca323e6b741
MIT
docs/contents/LJ_simulation.ipynb
uibcdf/OpenMembrane
c9705cb32706b882bdb3d75d19539ca323e6b741
MIT
Qwen/Qwen-72B
1. YES
2. YES