1#!/usr/local/bin/python3.8
2
3# DON'T RUN THIS MANUALLY
4# RUN "csg_inverse --options settings.xml" instead
5
6import time
7import numpy as np
8
9import espressomd
10from espressomd.io.writer import h5md
11
12def readgrofile(filename):
13    atompos=[]
14    box=[]
15    atomnumber=0
16    with open(filename,"r") as f:
17        f.readline()
18        atomnumber=int(f.readline())
19        for _ in range(atomnumber):
20            line=f.readline()
21            atompos.append(line.split()[3:])
22        box=(f.readline().split())
23    return atomnumber,np.array(box,dtype=float),np.array(atompos,dtype=float)
24
25def write_data(file_name, time, energy, n_part):
26    with open(file_name, 'a') as f:
27        e = energy['total'] / n_part
28        e_pot = (energy['total'] - energy['kinetic']) / n_part
29        e_kin = e - e_pot
30        print(" {} {} {} {}".format(time, e, e_pot, e_kin), file=f)
31
32
33def calc_temperature(system):
34    E = system.analysis.energy()['total']
35    n_part = len(system.part)
36
37    return 2. / 3. * E / n_part
38
39# check for necessary feature
40required_features = ["MASS", "TABULATED"]
41espressomd.assert_features(required_features)
42
43print("Program Information:\n{}\n".format(espressomd.features()))
44
45# set system properties:
46skin = 0.5
47time_step = 0.002
48friction = 5.0
49int_steps = 900
50eq_steps = 100
51steps_per_int = 100
52
53atomnumber,box,atompos=readgrofile('spce.gro')
54nm2Angstroem=10
55box_length = box*nm2Angstroem
56positions = atompos*nm2Angstroem
57mass=18
58masses=mass*np.ones(atomnumber)
59
60# Setup Espresso with particle from spce.gro
61system = espressomd.System(box_l=box_length, time_step=time_step)
62system.cell_system.skin = skin
63system.set_random_state_PRNG()
64system.thermostat.set_langevin(kT=1., gamma=1., seed=123)
65system.part.add(pos=positions, mass=masses)
66
67print("number of particles: ", len(system.part))
68print("box size: ", box_length)
69print("initial temperature: ", calc_temperature(system))
70
71with open('CG_CG.tab', 'r') as f:
72    data = np.loadtxt(f)
73    r = data[:, 0]
74    f = data[:, 1]
75    p = data[:, 2]
76
77    min_r = np.min(r)
78    max_r = np.max(r)
79
80    system.non_bonded_inter[0, 0].tabulated.set_params(min=min_r, max=max_r,
81                                                       energy=p,
82                                                       force=f)
83
84# Warmup loop
85for i in range(eq_steps):
86    print("warmup step", i)
87    system.integrator.run(steps_per_int)
88
89print("Running at temperature T={:.2f}".format(calc_temperature(system)))
90
91h5_file = h5md.H5md(filename="traj.h5", write_pos=True, write_vel=True,
92                    write_force=True, write_species=False, write_mass=True, write_ordered=False)
93
94starttime = time.time()
95
96# Integration loop
97for i in range(int_steps):
98    system.integrator.run(steps_per_int)
99
100    E = system.analysis.energy()
101    e_pot = (E['total'] - E['kinetic']) / len(system.part)
102    print("time: {:.3f} potential energy: {:.2f}".format(i * time_step, e_pot))
103
104    h5_file.write()
105    write_data("energy.dat", i * time_step,
106               system.analysis.energy(), len(system.part))
107
108endtime = time.time()
109h5_file.close()
110
111print("time per MD step: ", 1000. * (endtime - starttime) / int_steps, "ms")
112print("final temperature: ", calc_temperature(system))
113print("Finished.")
114