Dipole angle tutorial#
Here we analyse a box of water molecules exposed to an alternating electric field. To follow this tutorial, the data test files electricfwater of MAICoS are needed.
Simulation details - The system contains \(\sim 5000\) water molecules simulated for 100 ps in the NVT ensemble at a temperature of 300K. Periodic boundary conditions were employed in all directions and long range electrostatics were modelled using the PME method. LINCS algorithm was used to constraint the H-Bonds. The alternating electric field was applied along the \(x\)-axis in the form of a Gaussian laser pulse. Please check gromacs electric field for more details.
To follow this tutorial, the data test files airwater of MAICoS are needed. You can obtain them by cloning MAICoS repository:
git clone git@gitlab.com:maicos-devel/maicos.git
The electricfwater data files are located in tests/data/electricfwater/.
Option 1: from the Python interpreter#
First, let us ignore unnecessary warnings:
[1]:
import warnings
warnings.filterwarnings("ignore")
Then, import MAICoS, NumPy, MDAnalysis, and PyPlot:
[2]:
import maicos
import numpy as np
import MDAnalysis as mda
import matplotlib.pyplot as plt
from matplotlib.ticker import AutoMinorLocator
---------------------------------------------------------------------------
ImportError Traceback (most recent call last)
Cell In [2], line 1
----> 1 import maicos
2 import numpy as np
3 import MDAnalysis as mda
File ~/checkouts/readthedocs.org/user_builds/maicos/envs/v0.5/lib/python3.10/site-packages/maicos/__init__.py:36
29 from .modules.density import DensityPlanar, DensityCylinder
30 from .modules.epsilon import (
31 EpsilonBulk,
32 EpsilonPlanar,
33 EpsilonCylinder,
34 DielectricSpectrum
35 )
---> 36 from .modules.structure import Saxs, Debye, Diporder
37 from .modules.timeseries import DipoleAngle, KineticEnergy
38 from .modules.transport import Velocity
File ~/checkouts/readthedocs.org/user_builds/maicos/envs/v0.5/lib/python3.10/site-packages/maicos/modules/structure.py:24
18 from .. import tables
19 from ..decorators import (
20 make_whole,
21 set_planar_class_doc,
22 set_verbose_doc,
23 )
---> 24 from ..lib import sfactor
25 from ..utils import check_compound, savetxt
27 logger = logging.getLogger(__name__)
ImportError: dynamic module does not define module export function (PyInit_sfactor)
and set a few parameters for plotting purpose:
[3]:
fontsize = 25
font = {'family': 'sans', 'color': 'black',
'weight': 'normal', 'size': fontsize}
my_color_1 = np.array([0.090, 0.247, 0.560])
my_color_2 = np.array([0.235, 0.682, 0.639])
my_color_3 = np.array([1.000, 0.509, 0.333])
my_color_4 = np.array([0.588, 0.588, 0.588])
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [3], line 4
1 fontsize = 25
2 font = {'family': 'sans', 'color': 'black',
3 'weight': 'normal', 'size': fontsize}
----> 4 my_color_1 = np.array([0.090, 0.247, 0.560])
5 my_color_2 = np.array([0.235, 0.682, 0.639])
6 my_color_3 = np.array([1.000, 0.509, 0.333])
NameError: name 'np' is not defined
Define the path to the electricfwater data folder of MAICoS:
[4]:
datapath = "../../../../tests/data/electricfwater/"
Create a MDAnalysis universe, and extract its atoms:
[5]:
u = mda.Universe(datapath+"mdelectric.tpr",
datapath+"mdelectric.trr")
atoms = u.atoms
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [5], line 1
----> 1 u = mda.Universe(datapath+"mdelectric.tpr",
2 datapath+"mdelectric.trr")
3 atoms = u.atoms
NameError: name 'mda' is not defined
Important note: Because this calculation is computationaly heavy, let use isolate a small portion of the trajectory file. To perform the full calculation, just comment the following line:
[6]:
u.transfer_to_memory(stop = 100)
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [6], line 1
----> 1 u.transfer_to_memory(stop = 100)
NameError: name 'u' is not defined
Let us print a few information about the trajectory file:
[7]:
print(f"The number of water molecules is {np.int32(u.atoms.n_atoms/3)}")
timestep = np.round(u.trajectory.dt,2)
print(f"The time interval between the frames is {timestep} ps")
total_time = np.round(u.trajectory.totaltime,2)
print(f"The total simulation time is {total_time} ps")
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [7], line 1
----> 1 print(f"The number of water molecules is {np.int32(u.atoms.n_atoms/3)}")
2 timestep = np.round(u.trajectory.dt,2)
3 print(f"The time interval between the frames is {timestep} ps")
NameError: name 'np' is not defined
Run the MAICoS dipole angle function, selecting the \(x\) axis using the dim keyword:
[8]:
dipangle = maicos.DipoleAngle(atoms, dim=0)
dipangle.run()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [8], line 1
----> 1 dipangle = maicos.DipoleAngle(atoms, dim=0)
2 dipangle.run()
NameError: name 'maicos' is not defined
The run produces a python dictionary named dipangle.results with 4 keys linked to numpy arrays as values: 1. timestep, 2. cosine of dipole and x-axis, 3. cosine squared, 4. product of cosine of dipoles \(i\) and \(j\) with \(i \ne j\).
Extract the time vector and the \(\cos(\theta_i)\) data from the results attribute:
[9]:
t = dipangle.results["t"]
cos_theta_i = dipangle.results["cos_theta_i"]
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [9], line 1
----> 1 t = dipangle.results["t"]
2 cos_theta_i = dipangle.results["cos_theta_i"]
NameError: name 'dipangle' is not defined
Plot the final results using :
[10]:
fig = plt.figure(figsize=(13,6.5))
ax1 = plt.subplot(1, 1, 1)
plt.xlabel(r"$t$ (ps)", fontdict=font)
plt.ylabel(r"cos($\theta_i$)", fontdict=font)
plt.xticks(fontsize=fontsize)
plt.yticks(fontsize=fontsize)
ax1.plot(t, cos_theta_i, color=my_color_1, linewidth=4)
ax1.yaxis.offsetText.set_fontsize(20)
ax1.minorticks_on()
ax1.tick_params('both', length=10, width=2, which='major', direction='in')
ax1.tick_params('both', length=6, width=1.4, which='minor', direction='in')
ax1.xaxis.set_ticks_position('both')
ax1.yaxis.set_ticks_position('both')
ax1.spines["top"].set_linewidth(2)
ax1.spines["bottom"].set_linewidth(2)
ax1.spines["left"].set_linewidth(2)
ax1.spines["right"].set_linewidth(2)
ax1.yaxis.offsetText.set_fontsize(30)
minor_locator_y = AutoMinorLocator(2)
ax1.yaxis.set_minor_locator(minor_locator_y)
minor_locator_x = AutoMinorLocator(2)
ax1.xaxis.set_minor_locator(minor_locator_x)
ax1.tick_params(axis='x', pad=10)
plt.show()
---------------------------------------------------------------------------
NameError Traceback (most recent call last)
Cell In [10], line 1
----> 1 fig = plt.figure(figsize=(13,6.5))
2 ax1 = plt.subplot(1, 1, 1)
3 plt.xlabel(r"$t$ (ps)", fontdict=font)
NameError: name 'plt' is not defined
Option 2: from the command line interface#
Go to the electricfwater data folder of MAICoS, then use the maicos command directly from the terminal:
cd tests/data/electricfwater/
maicos DipoleAngle -s md.tpr -f md.trr -d 0
The output file dipangle.dat is similar to dipangle.results and contains the data in columns. To know the full list of options, have a look at the Inputs.