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