Commit c8382890 authored by Kai Szuttor's avatar Kai Szuttor

Merged espressomd/python.

parents 252c020a 46c5ec69
......@@ -51,3 +51,8 @@ doc/dg/background_errors.unsorted
# IDE
.idea
# Compiled python
*.pyc
.ycm_extra_conf.py
......@@ -231,7 +231,6 @@ set(CMAKE_CXX_FLAGS "${CMAKE_CXX_FLAGS} ${MPI_COMPILE_FLAGS}")
set(CMAKE_EXE_LINKER_FLAGS "${CMAKE_EXE_LINKER_FLAGS} ${MPI_LINK_FLAGS}")
include_directories(${MPI_INCLUDE_PATH})
list(APPEND LIBRARIES ${MPI_LIBRARIES})
add_definitions(-DH5XX_USE_MPI)
#######################################################################
# Paths
......
......@@ -15,5 +15,7 @@ add_custom_command(
# FIXME here: for the XX give it a proper number (01, 02...)
add_custom_target(tutorials_py_01 DEPENDS ${BASENAME}.pdf)
configure_file(scripts/lj_tutorial.py ${CMAKE_BINARY_DIR}/doc/tutorials/python/${BASENAME}/lj_tutorial.py COPYONLY)
configure_file(scripts/msd.py ${CMAKE_BINARY_DIR}/doc/tutorials/python/${BASENAME}/msd.py COPYONLY)
configure_file(scripts/lj_tutorial.py ${CMAKE_CURRENT_BINARY_DIR}/lj_tutorial.py COPYONLY)
configure_file(scripts/msd.py ${CMAKE_CURRENT_BINARY_DIR}/msd.py COPYONLY)
configure_file(scripts/two-component.py ${CMAKE_CURRENT_BINARY_DIR}/two-component.py COPYONLY)
configure_file(scripts/two-component-visualization.py ${CMAKE_CURRENT_BINARY_DIR}/two-component-visualization.py COPYONLY)
......@@ -20,7 +20,6 @@ from __future__ import print_function
import espressomd
from espressomd import code_info
import cPickle as pickle
import os
import numpy as np
......@@ -34,26 +33,26 @@ print(code_info.features())
# System parameters
#############################################################
n_part = 108
n_part = 500
density = 0.8442
skin = 0.1
time_step = 0.001
eq_tstep = 0.0001
time_step = 0.01
eq_tstep = 0.001
temperature = 0.728
box_l = np.power(n_part/density, 1.0/3.0) + 2*skin
box_l = np.power(n_part/density, 1.0/3.0)
warm_steps = 100
warm_n_time = 2000
min_dist = 0.87
# integration
sampling_interval = 10
sampling_interval = 100
equilibration_interval = 1000
sampling_iterations = 10000
equilibration_iterations= 20
sampling_iterations = 100
equilibration_iterations= 5
# Interaction parameters (Lennard Jones)
......@@ -62,7 +61,7 @@ equilibration_iterations= 20
lj_eps = 1.0
lj_sig = 1.0
lj_cut = 2.5*lj_sig
lj_cap = 20
lj_cap = 5
# System setup
......@@ -128,37 +127,24 @@ for i in range(equilibration_iterations):
print("eq run {} at time {}\n".format(i, system.time))
print("\nEquilibration done\n")
# Switch thermostat off NVE msd measurement
system.thermostat.turn_off()
print("\nSampling\n")
system.time_step = time_step
# Record energy versus time
en_fp = open('data/energy.dat', 'w')
sys_fp = open('data/sim_info.pickle', 'w')
part_fp = open('data/part.pickle', 'w')
msd_fp = open('data/msd.dat', 'w')
# Record radial distribution function
rdf_fp = open('data/rdf.dat', 'w')
en_fp.write("#\n#\n#\n# Pressure Kinetic Potential Temperature\n#\n")
#start positions of all particles for MSD calculation on the fly
start_pos=system.part[:].pos
en_fp.write("#\n#\n#\n# Time\ttotal energy\tkinetic energy\tlennard jones energy\ttemperature\n")
# save system setup
pickle.dump(system, sys_fp, -1)
msd = np.zeros((sampling_iterations,))
# Data arrays for simple error estimation
etotal = np.zeros((sampling_iterations,))
ptotal = np.zeros((sampling_iterations,))
# save start particle configuration
pickle.dump(system.part, part_fp, -1)
# analyzing the radial distribution function
# setting the parameters for the rdf
r_bins = 30
r_bins = 50
r_min = 0.0
r_max = system.box_l[0]/2.0
......@@ -167,45 +153,27 @@ avg_rdf=np.zeros((r_bins,))
for i in range(1, sampling_iterations + 1):
system.integrator.run(sampling_interval)
energies = system.analysis.energy()
pressure = system.analysis.pressure()
r, rdf = system.analysis.rdf(rdf_type="rdf", type_list_a=[0], type_list_b=[0], r_min=r_min, r_max=r_max, r_bins=r_bins)
avg_rdf+= rdf/sampling_iterations
kinetic_temperature = energies['ideal']/( 1.5 * n_part)
en_fp.write("%i\t%1.5e\t%1.5e\t%1.5e\t%1.5e\t%1.5e\n" % (i, pressure['total'], energies['total'], energies['ideal'], energies['total'] - energies['ideal'], kinetic_temperature))
en_fp.write("%f\t%1.5e\t%1.5e\t%1.5e\t%1.5e\n" % (system.time, energies['total'], energies['ideal'], energies['total'] - energies['ideal'], kinetic_temperature))
etotal[i-1] = energies['total']
ptotal[i-1] = pressure['total']
####################################################
# TODO:
# For an efficient calculation of the MSD and VACF the correlator and
# observable features are needed in python
####################################################
for j in range(n_part):
msd[i-1] += 1.0/n_part * ( (start_pos[j,0] - system.part[j].pos[0])**2 +
(start_pos[j,1] - system.part[j].pos[1])**2 +
(start_pos[j,2] - system.part[j].pos[2])**2 )
msd_fp.write("%i %1.5e\n" % (i, msd[i-1]) )
print("\nMain sampling done\n")
# calculate the variance of the total energy and total pressure using scipys statistic operations
# calculate the variance of the total energy using scipys statistic operations
error_total_energy=np.sqrt(etotal.var())/np.sqrt(sampling_iterations)
error_total_pressure=np.sqrt(ptotal.var())/np.sqrt(sampling_iterations)
en_fp.write("#mean_energy energy_error mean_pressure pressure_error\n#%1.5e %1.5e %1.5e %1.5e" \
% (etotal.mean(), error_total_energy, ptotal.mean(), error_total_pressure) )
en_fp.write("#mean_energy energy_error %1.5e %1.5e\n" % (etotal.mean(), error_total_energy) )
# write out the radial distribution data
for i in range(r_bins):
rdf_fp.write("%1.5e %1.5e\n" % (r[i], avg_rdf[i]))
sys_fp.close()
msd_fp.close()
en_fp.close()
part_fp.close()
rdf_fp.close()
......@@ -16,37 +16,183 @@
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# This is a modified version of the lj liquid script, which calculates the
# mean saure displacement versus time.
# 1. Setup up and equilibrate the Lennard Jones liquid
# 2. Setup an observable for the particle positions
# 3. Set up a "Correlator" which records the square distance of the particles
# for different time intervals
# 4. Integrate equations of motion
# 5. Print the result from the correlator
# 1. Setup and equilibrate LJ liquid
from __future__ import print_function
import espressomd._system as es
import espressomd
from espressomd import thermostat
from espressomd import code_info
from espressomd import integrate
import cPickle as pickle
import os
import numpy as np
data_file="data/sim_info.pickle"
part_file="data/part.pickle"
n_part = 200
density = 0.8442
skin = 0.1
time_step = 0.01
eq_tstep = 0.01
temperature = 0.728
box_l = np.power(n_part/density, 1.0/3.0)
warm_steps = 100
warm_n_time = 2000
min_dist = 0.87
# integration
sampling_interval = 10
equilibration_interval = 1000
sampling_iterations = 10000
equilibration_iterations= 10
# Interaction parameters (Lennard Jones)
#############################################################
lj_eps = 1.0
lj_sig = 1.0
lj_cut = 2.5*lj_sig
lj_cap = 20
# System setup
#############################################################
system = espressomd.System()
if not os.path.exists('data') :
os.mkdir('data')
system.time_step = time_step
system.cell_system.skin = skin
system.box_l = [box_l, box_l, box_l]
system.non_bonded_inter[0, 0].lennard_jones.set_params(
epsilon=lj_eps, sigma=lj_sig,
cutoff=lj_cut, shift="auto")
system.non_bonded_inter.set_force_cap(lj_cap)
print("LJ-parameters:")
print(system.non_bonded_inter[0, 0].lennard_jones.get_params())
# Thermostat
system.thermostat.set_langevin(kT=temperature, gamma=1.0)
# Particle setup
#############################################################
volume = box_l * box_l * box_l
system = espressomd.System()
pickle.load(open(data_file,"r"))
for i in range(n_part):
system.part.add(id=i, pos=np.random.random(3) * system.box_l)
#############################################################
# Warmup Integration #
#############################################################
print("""
ro variables:
cell_grid {0.cell_grid}
cell_size {0.cell_size}
local_box_l {0.local_box_l}
max_cut {0.max_cut}
max_part {0.max_part}
max_range {0.max_range}
max_skin {0.max_skin}
n_nodes {0.n_nodes}
n_part {0.n_part}
n_part_types {0.n_part_types}
periodicity {0.periodicity}
transfer_rate {0.transfer_rate}
""".format(system))
pickle.load(open(part_file,"r"))
print(system.part)
Start warmup integration:
At maximum {} times {} steps
Stop if minimal distance is larger than {}
""".strip().format(warm_n_time, warm_steps, min_dist))
i = 0
act_min_dist = system.analysis.mindist()
while i < warm_n_time and act_min_dist < min_dist :
system.integrator.run(warm_steps)
act_min_dist = system.analysis.mindist()
print("run {} at time = {} (LJ cap= {} ) min dist = {}".strip().format(i, system.time, lj_cap, act_min_dist))
i+=1
lj_cap += 1.0
system.non_bonded_inter.set_force_cap(lj_cap)
system.non_bonded_inter.set_force_cap(0)
print("\nWarm up finished\n")
system.time_step = eq_tstep
for i in range(equilibration_iterations):
system.integrator.run(equilibration_interval)
energies = system.analysis.energy()
print("eq run {} at time {}\n".format(i, system.time))
print("\nEquilibration done\n")
system.time_step = time_step
# 2. Setup an observable for the particle positions
from espressomd.observables import ParticlePositions
# We pass the ids of the particles to be tracked to the observable.
# Here this is 0..n_part
part_pos=ParticlePositions(ids=range(n_part))
# 3. Define a correlator recording the square distances the particles have traveled
# in various time intervals
from espressomd.correlators import Correlator
# The correlator works with the part_pos observable. Here, no second
# observableis needed
# For short time spans, we record in linear time distances
# for larger time spans, we record in larger and larger steps.
# Use 10 linear measurements 10 time steps apart, each.
# The "square_distance_component_wise" correlation operation tells
# the correlator how to calculate a result from the measurements of the
# observables at different times
corr=Correlator(obs1=part_pos,
tau_lin=10,dt=10*time_step,
tau_max=10000*time_step,
corr_operation="square_distance_componentwise")
# We want the correlator to take measurements and calculate results
# automatically during the integration
system.auto_update_correlators.add(corr)
# 4. Integrate the equations of motion
# Note that a longer integration would be needed to achieve good quality
# results. This will run for ~5 minutes.
system.integrator.run(100000)
# 5. Save the result of the correlation to a file
# The 1st column contains the time, the 2nd column the number of
# measurements, and the remaining columns the mean square displacement
# for each particle and each Cartesian component
corr.finalize()
result=corr.result()
# We average over all particls and coordinaes and store id in the 3rd column
# using numpy's averaging function
result[:,3]=np.average(result[:,3:],1)
# And save the first three columns to a file
np.savetxt("data/msd.dat",result[:,:3])
# Plot the msd (third against first column of msd.dat) with your favorite
# plotting program and measure the diffusion coefficient by fitting
# <x^2> =2Dt
#
# Copyright (C) 2013,2014,2015,2016 The ESPResSo project
#
# This file is part of ESPResSo.
#
# ESPResSo 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 3 of the License, or
# (at your option) any later version.
#
# ESPResSo 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.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# This is a modified version of the lj liquid script, which simulates
# a two component lj liquid.
# By switching the lj interaction between the two components
# from atractive to purely repulsive, de-mixing can be achieved.
# 1. Setup and equilibrate LJ liquid
from __future__ import print_function
import espressomd
from espressomd import code_info
import os
import numpy as np
from espressomd import visualization
from threading import Thread
n_part = 200
density = 0.4442
skin = 0.1
time_step = 0.01
eq_tstep = 0.01
temperature = 0.728
box_l = np.power(n_part/density, 1.0/3.0)
warm_steps = 100
warm_n_time = 2000
min_dist = 0.87
# integration
sampling_interval = 10
equilibration_interval = 1000
sampling_iterations = 10000
equilibration_iterations= 10
# Interaction parameters (Lennard Jones)
#############################################################
lj_eps = 1.0
lj_sig = 1.0
lj_cut = 2.5*lj_sig
lj_cap = 20
# This is the cutoff of the interaction between species 0 and 1.
# By setting it to 2**(1./6.) *lj_sig, it can be made purely repulsive
lj_cut_mixed =2.5 * lj_sig
lj_cut_mixed =2**(1./6.) * lj_sig
# System setup
#############################################################
system = espressomd.System()
if not os.path.exists('data') :
os.mkdir('data')
system.time_step = time_step
system.cell_system.skin = skin
system.box_l = [box_l, box_l, box_l]
# Here, lj interactions need to be setup for both components
# as well as for the mixed case of component 0 interacting with
# component 1.
# component 0
system.non_bonded_inter[0, 0].lennard_jones.set_params(
epsilon=lj_eps, sigma=lj_sig,
cutoff=lj_cut, shift="auto")
# component 1
system.non_bonded_inter[1, 1].lennard_jones.set_params(
epsilon=lj_eps, sigma=lj_sig,
cutoff=lj_cut, shift="auto")
# mixed case
system.non_bonded_inter[0, 1].lennard_jones.set_params(
epsilon=lj_eps, sigma=lj_sig,
cutoff=lj_cut_mixed, shift="auto")
system.non_bonded_inter.set_force_cap(lj_cap)
print("LJ-parameters:")
print(system.non_bonded_inter[0, 0].lennard_jones.get_params())
# Thermostat
system.thermostat.set_langevin(kT=temperature, gamma=1.0)
# Particle setup
#############################################################
volume = box_l * box_l * box_l
for i in range(n_part):
system.part.add(id=i, pos=np.random.random(3) * system.box_l)
# Every 2nd particle should be of component 1
if i%2==1: system.part[i].type=1
#############################################################
# Warmup Integration #
#############################################################
print("""
Start warmup integration:
At maximum {} times {} steps
Stop if minimal distance is larger than {}
""".strip().format(warm_n_time, warm_steps, min_dist))
i = 0
act_min_dist = system.analysis.mindist()
while i < warm_n_time and act_min_dist < min_dist :
system.integrator.run(warm_steps)
act_min_dist = system.analysis.mindist()
print("run {} at time = {} (LJ cap= {} ) min dist = {}".strip().format(i, system.time, lj_cap, act_min_dist))
i+=1
lj_cap += 1.0
system.non_bonded_inter.set_force_cap(lj_cap)
system.non_bonded_inter.set_force_cap(0)
def loop():
while True:
system.integrator.run(100)
visualizer.update()
visualizer = visualization.mayaviLive(system)
#Start simulation in seperate thread
t = Thread(target=loop)
t.daemon = True
t.start()
#Start blocking visualizer
visualizer.start()
#
# Copyright (C) 2013,2014,2015,2016 The ESPResSo project
#
# This file is part of ESPResSo.
#
# ESPResSo 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 3 of the License, or
# (at your option) any later version.
#
# ESPResSo 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.
#
# You should have received a copy of the GNU General Public License
# along with this program. If not, see <http://www.gnu.org/licenses/>.
#
# This is a modified version of the lj liquid script, which simulates
# a two component lj liquid.
# By switching the lj interaction between the two components
# from atractive to purely repulsive, de-mixing can be achieved.
# 1. Setup and equilibrate LJ liquid
from __future__ import print_function
import espressomd
from espressomd import code_info
import os
import numpy as np
n_part = 200
density = 0.8442
skin = 0.1
time_step = 0.01
eq_tstep = 0.01
temperature = 0.728
box_l = np.power(n_part/density, 1.0/3.0)
warm_steps = 100
warm_n_time = 2000
min_dist = 0.87
# integration
sampling_interval = 10
equilibration_interval = 1000
sampling_iterations = 10000
equilibration_iterations= 10
# Interaction parameters (Lennard Jones)
#############################################################
lj_eps = 1.0
lj_sig = 1.0
lj_cut = 2.5*lj_sig
lj_cap = 20
# This is the cutoff of the interaction between species 0 and 1.
# By setting it to 2**(1./6.) *lj_sig, it can be made purely repulsive
lj_cut_mixed =2.5 * lj_sig
lj_cut_mixed =2**(1./6.) * lj_sig
# System setup
#############################################################
system = espressomd.System()
if not os.path.exists('data') :
os.mkdir('data')
system.time_step = time_step
system.cell_system.skin = skin
system.box_l = [box_l, box_l, box_l]
# Here, lj interactions need to be setup for both components
# as well as for the mixed case of component 0 interacting with
# component 1.
# component 0
system.non_bonded_inter[0, 0].lennard_jones.set_params(
epsilon=lj_eps, sigma=lj_sig,
cutoff=lj_cut, shift="auto")
# component 1
system.non_bonded_inter[1, 1].lennard_jones.set_params(
epsilon=lj_eps, sigma=lj_sig,
cutoff=lj_cut, shift="auto")
# mixed case
system.non_bonded_inter[0, 1].lennard_jones.set_params(
epsilon=lj_eps, sigma=lj_sig,
cutoff=lj_cut_mixed, shift="auto")
system.non_bonded_inter.set_force_cap(lj_cap)
print("LJ-parameters:")
print(system.non_bonded_inter[0, 0].lennard_jones.get_params())
# Thermostat
system.thermostat.set_langevin(kT=temperature, gamma=1.0)
# Particle setup
#############################################################
volume = box_l * box_l * box_l
for i in range(n_part):
system.part.add(id=i, pos=np.random.random(3) * system.box_l)
# Every 2nd particle should be of component 1
if i%2==1: system.part[i].type=1
#############################################################
# Warmup Integration #
#############################################################
print("""
Start warmup integration:
At maximum {} times {} steps
Stop if minimal distance is larger than {}
""".strip().format(warm_n_time, warm_steps, min_dist))
i = 0
act_min_dist = system.analysis.mindist()
while i < warm_n_time and act_min_dist < min_dist :
system.integrator.run(warm_steps)
act_min_dist = system.analysis.mindist()
print("run {} at time = {} (LJ cap= {} ) min dist = {}".strip().format(i, system.time, lj_cap, act_min_dist))
i+=1
lj_cap += 1.0
system.non_bonded_inter.set_force_cap(lj_cap)
system.non_bonded_inter.set_force_cap(0)
print("\nWarm up finished\n")
system.time_step = eq_tstep
for i in range(equilibration_iterations):
system.integrator.run(equilibration_interval)
energies = system.analysis.energy()