espressomd-users
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

Re: [ESPResSo-users] writing log and trajectories


From: Kai Szuttor
Subject: Re: [ESPResSo-users] writing log and trajectories
Date: Fri, 6 Sep 2019 13:15:11 +0200

Hi,

could you please tell us which version of ESPReSso you are using?
If you are using the current release, find documentation for file I/O here: http://espressomd.org/html/doc4.0.2/io.html

Best,

Kai


Kai Szuttor
Institute for Computational Physics
Universität Stuttgart
Allmandring 3
Room 1.084
70569 Stuttgart, Germany
Phone: +49 711 685-67707

On 6. Sep 2019, at 13:03, Raghvendra Singh <address@hidden> wrote:

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')
#########################################################################################################################################


reply via email to

[Prev in Thread] Current Thread [Next in Thread]