[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[ESPResSo-users] writing log and trajectories
From: |
Raghvendra Singh |
Subject: |
[ESPResSo-users] writing log and trajectories |
Date: |
Fri, 6 Sep 2019 13:03:31 +0200 (CEST) |
Hello All,
I am newbie to EspressoMD and in a testing trial phase of it.
I am trying to run a simple LJ simulation with multispecies simulation system.
There is a problem with writing appropriate log and trajectory for my
simulation.
I went through the mailing list archive and documentation and still find it
hard to manage what I want.
this is a script I am using to run the simulation, please help me out with
getting the appropriate outputs required for further analysis.
I would to have an output log and a trajectory every 100 step.
Thank you all
Best
Raghav
#########################################################################################################################################
import espressomd
import MDAnalysis as mda
from espressomd.visualization_opengl import openGLLive
from espressomd import electrostatics
from espressomd.io.writer import vtf
from espressomd import MDA_ESP
import numpy as np
required_features = ["P3M", "LENNARD_JONES", "MASS"]
espressomd.assert_features(required_features)
box = [200, 200, 200]
system = espressomd.System(box_l=box)
fp=open('trajectory.vtf', mode='w+t')
system.cell_system.set_domain_decomposition(use_verlet_lists=True)
#visualizer = openGLLive(system, background_color=[1, 1, 1],
drag_enabled=True, drag_force=10)
system.set_random_state_PRNG()
system.seed=42
# TIMESTEP
time_step_fs = 1.0
system.time_step = time_step_fs * 1.0e-2
system.cell_system.skin = 1.2
# TEMPERATURE
SI_temperature = 300.0
kb_kjmol = 0.0083145
temperature = SI_temperature * kb_kjmol
# COULOMB PREFACTOR (elementary charge)^2 / (4*pi*epsilon_0) in Angstrom *
# kJ/mol
epsilon_r = 4.0
coulomb_prefactor = 1.67101e5 * kb_kjmol / epsilon_r
# FORCE FIELDS
species = ["Cl", "Na", "C1", "C2","C3","C4","CCP","Solvent"]
types = {"Cl": 0, "Na": 1, "C1": 2,"C2":3,"C3":4,"C4":5,"CCP":6,"Solvent":7}
charges = {"Cl": -1.0, "Na": 1.0, "C1": -6.0, "C2": 0.0, "C3": 0.0, "C4":
0.0,"CCP":6,"Solvent": 0.0}
lj_sigmas = {"Cl": 3.85, "Na": 2.52, "C1": 10.0, "C2": 10.0, "C3": 10.0, "C4":
10.0,"CCP":5, "Solvent": 1.5}
lj_epsilons = {"Cl": 192.45, "Na": 17.44,
"C1": 100.0, "C2": 100.0,
"C3": 100.0, "C4": 100.0,
"CCP": 70, "Solvent": 50.0}
lj_cuts = {"Cl": 2.0 * lj_sigmas["Cl"], "Na": 2.0 * lj_sigmas["Na"],
"C1": 5.0 * lj_sigmas["C1"],"C2": 5.0 * lj_sigmas["C2"],
"C3": 5.0 * lj_sigmas["C3"],"C4": 5.0 * lj_sigmas["C4"],
"CCP": 2.5 * lj_sigmas["CCP"],"Solvent": 2.0 * lj_sigmas["Solvent"]}
masses = {"Cl": 35.453, "Na": 22.99, "C1": 100, "C2": 100, "C3": 100, "C4":
100,"CCP": 50, "Solvent": 18.0}
n_ionpairs = 150
for i in range(n_ionpairs):
for t in ["Na", "Cl"]:
system.part.add(pos=box * np.random.random(3),
q=charges[t], type=types[t], mass=masses[t])
n_colloids = 150
t = "C1"
t_co = "Cl"
for i in range(n_colloids):
system.part.add(pos=box * np.random.random(3),
q=charges[t], type=types[t], mass=masses[t])
for i in range(int(abs(charges[t]))):
system.part.add(pos=box * np.random.random(3),
q=charges[t_co], type=types[t_co], mass=masses[t_co])
n_colloids = 150
t = "C2"
t_co = "Na"
for i in range(n_colloids):
system.part.add(pos=box * np.random.random(3),
q=charges[t], type=types[t], mass=masses[t])
for i in range(int(abs(charges[t]))):
system.part.add(pos=box * np.random.random(3),
q=charges[t_co], type=types[t_co], mass=masses[t_co])
n_colloids = 150
t = "C3"
t_co = "Cl"
for i in range(n_colloids):
system.part.add(pos=box * np.random.random(3),
q=charges[t], type=types[t], mass=masses[t])
for i in range(int(abs(charges[t]))):
system.part.add(pos=box * np.random.random(3),
q=charges[t_co], type=types[t_co], mass=masses[t_co])
n_colloids = 150
t = "C4"
t_co = "Cl"
for i in range(n_colloids):
system.part.add(pos=box * np.random.random(3),
q=charges[t], type=types[t], mass=masses[t])
for i in range(int(abs(charges[t]))):
system.part.add(pos=box * np.random.random(3),
q=charges[t_co], type=types[t_co], mass=masses[t_co])
n_colloids = 150
t = "CCP"
t_co = "Na"
for i in range(n_colloids):
system.part.add(pos=box * np.random.random(3),
q=charges[t], type=types[t], mass=masses[t])
for i in range(int(abs(charges[t]))):
system.part.add(pos=box * np.random.random(3),
q=charges[t_co], type=types[t_co], mass=masses[t_co])
n_solvents = 0
t = "Solvent"
for i in range(n_solvents):
system.part.add(pos=box * np.random.random(3),
q=charges[t], type=types[t], mass=masses[t])
def combination_rule_epsilon(rule, eps1, eps2):
if rule == "Lorentz":
return (eps1 * eps2)**0.5
else:
return ValueError("No combination rule defined")
def combination_rule_sigma(rule, sig1, sig2):
if rule == "Berthelot":
return (sig1 + sig2) * 0.5
else:
return ValueError("No combination rule defined")
# Lennard-Jones interactions parameters
for i in range(len(species)):
for j in range(i, len(species)):
s = [species[i], species[j]]
lj_sig = combination_rule_sigma(
"Berthelot", lj_sigmas[s[0]], lj_sigmas[s[1]])
lj_cut = combination_rule_sigma(
"Berthelot", lj_cuts[s[0]], lj_cuts[s[1]])
lj_eps = combination_rule_epsilon(
"Lorentz", lj_epsilons[s[0]], lj_epsilons[s[1]])
system.non_bonded_inter[types[s[0]],
types[s[1]]].lennard_jones.set_params(
epsilon=lj_eps, sigma=lj_sig, cutoff=lj_cut, shift="auto")
energy = system.analysis.energy()
print("Before Minimization: E_total = {}".format(energy['total']))
system.minimize_energy.init(f_max=1000, gamma=30.0,
max_steps=10000, max_displacement=0.01)
system.minimize_energy.minimize()
energy = system.analysis.energy()
print("After Minimization: E_total = {}".format(energy['total']))
print("Tune p3m")
p3m = electrostatics.P3M(prefactor=coulomb_prefactor, accuracy=1e-1)
system.actors.add(p3m)
system.thermostat.set_langevin(kT=temperature, gamma=2.0)
vtf.writevsf(system, fp, types='all')
vtf.writevcf(system, fp, types='all')
for n in num_steps:
system.integrator.run(100)
vtf.writevcf(system, fp)
fp.close()
vtf.writevsf(system, fp, types='all')
#########################################################################################################################################
- [ESPResSo-users] writing log and trajectories,
Raghvendra Singh <=