import espressomd
import espressomd.observables
import espressomd.accumulators
import espressomd.analyze
import espressomd.interactions
import numpy as np
required_features = ["DPD", "LENNARD_JONES", "LJCOS2"]
espressomd.assert_features(required_features)
mol_num = 320
box_l = [14.0, 14.0, 14.0]
system = espressomd.System(box_l=box_l)
# Set bonded interactions
fene = espressomd.interactions.FeneBond(k=30, d_r_max=1.5)
hb = espressomd.interactions.HarmonicBond(k=10.0, r_0=4.0)
# bonded_parms = [fene, hb]
system.bonded_inter.add(fene)
system.bonded_inter.add(hb)
# Set non-bonded interactions
lj_eps = 1.0
lj_sigmah = 0.95
lj_sigma = 1.0
lj_offset = 0.0
system.non_bonded_inter[0, 0].lennard_jones.set_params(
epsilon=lj_eps, sigma=lj_sigmah, cutoff=1.1225 * lj_sigmah,
shift=0.25 * lj_eps, offset=lj_offset)
system.non_bonded_inter[0, 1].lennard_jones.set_params(
epsilon=lj_eps, sigma=lj_sigmah, cutoff=1.1225 * lj_sigmah,
shift=0.25 * lj_eps, offset=lj_offset)
system.non_bonded_inter[1, 1].lennard_jones_cos2.set_params(
epsilon=lj_eps, sigma=lj_sigma, offset=lj_offset, width=1.6)
particle_types = [0, 1, 1]
bond_l = 1.0
for i in range(mol_num):
tail_pos = np.random.random(3) * system.box_l
orient = 2 * np.random.random(3) - 1
orient /= np.linalg.norm(orient)
for j in range(len(particle_types)):
cur_part_id = i * len(particle_types) + j
particle_position = tail_pos + (len(particle_types) - j - 1) * bond_l * orient
system.part.add(id=cur_part_id, type=particle_types[j], pos=particle_position, rotation=(True, True, True))
if j > 0:
system.part.by_id(cur_part_id - 1).add_bond((fene, cur_part_id))
if j > 1:
system.part.by_id(cur_part_id - 2).add_bond((hb, cur_part_id))
# Integration parameters
langevin_gamma = 1.0 # not inited
main_time_step = 0.01
verlet_skin = 0.4
systemtemp = 1.1
dpd_gamma = 1.0
dpd_r_cut = 1.12 + 1.8
int_steps = 200
int_n_times = 1000
def warmup(steps, times, startcap = 5, capgoal = 1000, min_dist = 0):
capincr = (capgoal - startcap) / times
cap = startcap
for i in range(times):
act_min_dist = system.analysis.min_dist()
if act_min_dist < min_dist:
break
system.force_cap = cap
system.integrator.run(steps)
cap += capincr
system.force_cap = 0
# Warmup parameters
warm_time_step = 0.002
warmup_temp = 0
warmsteps = 100
warmtimes = 20
free_warmsteps = 0
free_warmtimes = 1
# Warmup
if warmup_temp == 0:
warmup_temp = systemtemp
system.time_step = warm_time_step
system.cell_system.skin = verlet_skin
system.thermostat.set_langevin(kT=warmup_temp, gamma=langevin_gamma, seed=42)
warmup(warmsteps, warmtimes)
# Second warmup
system.time_step = main_time_step
system.thermostat.set_langevin(kT=systemtemp, gamma=langevin_gamma, seed=42)
warmup(free_warmsteps, free_warmtimes, 1000)
system.time = 0
# Integration
system.thermostat.turn_off()
system.thermostat.set_dpd(kT=systemtemp, seed=41)
system.non_bonded_inter[0, 0].dpd.set_params(gamma=dpd_gamma, r_cut=dpd_r_cut)
system.non_bonded_inter[0, 1].dpd.set_params(gamma=dpd_gamma, r_cut=dpd_r_cut)
system.non_bonded_inter[1, 1].dpd.set_params(gamma=dpd_gamma, r_cut=dpd_r_cut)
for i in range(int_n_times):
print("\rrun %d at time=%.0f " % (i, system.time), end='')
system.integrator.run(int_steps)