This code requires the gsd
package, which is available on conda forge.
Coarse-graining a box of 1000 methanol molecules from the all atom simulation given in CG_tutorial/meth.gsd
file.
# disable GPU. Remove this if you've compiled HOOMD for GPU
import os
os.environ['CUDA_VISIBLE_DEVICES'] = '-1'
import tensorflow as tf
import hoomd
import hoomd.md
import hoomd.htf as htf
import numpy as np
import gsd, gsd.hoomd
import matplotlib.pyplot as plt,matplotlib
Here we prepare the computations that will be executed at each step during the simulation. We have access to the neighbor list, positions, types, box dimensions of the simulation.
We work in four steps:
# set-up the system
g = gsd.hoomd.open('CG_tutorial/meth.gsd')
c = hoomd.context.initialize('--mode=cpu')
system = hoomd.init.read_gsd(filename='CG_tutorial/meth.gsd')
# disabled particle sorting!
c.sorter.disable()
set_rcut = 10.0
HOOMD-blue 2.5.2 DOUBLE HPMC_MIXED TBB SSE SSE2 SSE3 Compiled: 04/30/2019 Copyright (c) 2009-2019 The Regents of the University of Michigan. ----- You are using HOOMD-blue. Please cite the following: * J A Anderson, C D Lorenz, and A Travesset. "General purpose molecular dynamics simulations fully implemented on graphics processing units", Journal of Computational Physics 227 (2008) 5342--5359 * J Glaser, T D Nguyen, J A Anderson, P Liu, F Spiga, J A Millan, D C Morse, and S C Glotzer. "Strong scaling of general-purpose molecular dynamics simulations on GPUs", Computer Physics Communications 192 (2015) 97--107 ----- HOOMD-blue is running on the CPU notice(2): Group "all" created containing 6000 particles
#get mapping from molecule index to particle index
molecule_mapping_index = htf.find_molecules(system)
# get number of atoms
N = sum([len(m) for m in molecule_mapping_index])
# get number of molecules
M = len(molecule_mapping_index)
# get number of atoms in a molecule=MN
MN = len(molecule_mapping_index[0])
## TRY CHANGING DIFFERENT BEADS DISTRIBUTION FOR TESTING
# create one bead mapping ->
# 1 x 6: [1, 1, 1, 1, 1, 1] that means
# all atoms contribute to CG bead equally
# massess are accounted for in sparse_mapping
# molecule_mapping = np.ones([1, MN], dtype=np.int)
# create 2 bead mapping ->
# 2 x 6:
# [0, 0, 1, 1, 1, 1]
# [1, 1, 0, 0, 0, 0]
# 4 atoms in first bead, 2 in second
molecule_mapping = np.array([[1, 1, 1, 0, 0, 0], [0, 0, 0, 1, 1, 1]])
bead_number = molecule_mapping.shape[0]
#create a mass-weighted M x N mapping operator
cg_mapping = htf.sparse_mapping([molecule_mapping for _ in molecule_mapping_index],
molecule_mapping_index, system=system)
assert cg_mapping.shape == (M * bead_number, N)
Finding molecules...99.90%
# define model
class MappingModel(htf.SimModel):
def setup(self, CG_NN, cg_mapping, rcut):
self.CG_NN = CG_NN
self.rcut = rcut
self.cg_mapping = cg_mapping
self.avg_cg_rdf = tf.keras.metrics.MeanTensor()
self.avg_aa_rdf = tf.keras.metrics.MeanTensor()
def compute(self, nlist, positions, box):
# calculate the center of mass of a CG bead
box_size = htf.box_size(box)
mapped_pos = htf.center_of_mass(positions[:,:3], self.cg_mapping, box_size)
# create the mapped neighbot list
mapped_nlist = htf.compute_nlist(mapped_pos, self.rcut, self.CG_NN, box_size, True)
# compute RDF for mapped and C-C in all-atom
cg_rdf = htf.compute_rdf(mapped_nlist, [0.1,self.rcut])
aa_rdf = htf.compute_rdf(nlist, [0.1,self.rcut], positions[:,3], type_i=3, type_j=3)
self.avg_cg_rdf.update_state(cg_rdf)
self.avg_aa_rdf.update_state(aa_rdf)
return
# build model. output_forces = False because
# this model doesn't comptue forces
max_neighbor_est = 256
model = MappingModel(
max_neighbor_est,
CG_NN=max_neighbor_est,
cg_mapping=cg_mapping,
output_forces=False,
rcut=set_rcut,
check_nlist=True)
Here we begin a simulation of methanol. This code is a little complex, but not really about hoomd-tf. This is the details of setting-up the force field to simulate methanol (e.g., treating electrostatics, dispersion, thermostat).
### Hoomd-Sim code ###
tfcompute = htf.tfcompute(model)
nlist = hoomd.md.nlist.cell()
#set-up pppm
charged = hoomd.group.all()
pppm = hoomd.md.charge.pppm(nlist=nlist, group=charged)
pppm.set_params(Nx=32, Ny=32, Nz=32, order=6, rcut=set_rcut)
#set-up pair coefficients
nlist.reset_exclusions(['1-2', '1-3', '1-4','body'])
lj = hoomd.md.pair.force_shifted_lj(r_cut=set_rcut, nlist=nlist)
forces = [lj]
lj.pair_coeff.set("opls_156", "opls_156", sigma=2.5, epsilon=0.03)
lj.pair_coeff.set("opls_156", "opls_157", sigma=2.96, epsilon=0.05)
lj.pair_coeff.set("opls_156", "opls_154", sigma=2.79, epsilon=0.07)
lj.pair_coeff.set("opls_156", "opls_155", sigma=5.0, epsilon=0.0)
lj.pair_coeff.set("opls_157", "opls_157", sigma=3.5, epsilon=0.07)
lj.pair_coeff.set("opls_157", "opls_154", sigma=3.31, epsilon=0.11)
lj.pair_coeff.set("opls_157", "opls_155", sigma=5.92, epsilon=0.0)
lj.pair_coeff.set("opls_154", "opls_154", sigma=3.12, epsilon=0.17)
lj.pair_coeff.set("opls_154", "opls_155", sigma=5.59, epsilon=0.0)
lj.pair_coeff.set("opls_155", "opls_155", sigma=10.0, epsilon=0.0)
#set-up bonds
harmonic = hoomd.md.bond.harmonic()
harmonic.bond_coeff.set("opls_156-opls_157", k=340.00, r0=1.09)
harmonic.bond_coeff.set("opls_154-opls_157", k=320.00, r0=1.41)
harmonic.bond_coeff.set("opls_154-opls_155", k=553.00, r0=0.95)
#set-up angles
harm_angle = hoomd.md.angle.harmonic()
harm_angle.angle_coeff.set("opls_154-opls_157-opls_156", k=70.0, t0=1.90)
harm_angle.angle_coeff.set("opls_155-opls_154-opls_157", k=110.0, t0=1.89)
harm_angle.angle_coeff.set("opls_156-opls_157-opls_156", k=66.0, t0=1.88)
#set-up dihedrals
dihedral = hoomd.md.dihedral.opls()
dihedral.dihedral_coeff.set("opls_155-opls_154-opls_157-opls_156", k1=0.0, k2=0.0, k3=0.45, k4=0.0)
group_all = hoomd.group.all()
kT = 1.9872/1000
#NVT Simulation in Hoomd
im = hoomd.md.integrate.mode_standard(dt=5.0/489.0)
nvt = hoomd.md.integrate.nvt(group=group_all, kT=298.15 * kT, tau=350 / 48.9)
nvt.randomize_velocities(1234)
tfcompute.attach(nlist, r_cut=set_rcut)
#Hoomd production run
hoomd.run(250)
----- You are using PPPM. Please cite the following: * D N LeBard, B G Levine, S A Barr, A Jusufi, S Sanders, M L Klein, and A Z Panagiotopoulos. "Self-assembly of coarse-grained ionic surfactants accelerated by graphics processing units", Journal of Computational Physics 8 (2012) 2385-2397 ----- notice(2): -- Neighborlist exclusion statistics -- : notice(2): Particles with 5 exclusions : 6000 notice(2): Neighbors included by diameter : no notice(2): Neighbors excluded when in the same body: no notice(2): Force mode is FORCE_MODE.hoomd2tf notice(2): Starting TensorflowCompute notice(2): completed reallocate ** starting run **
*Warning*: charge.pppm: system is not neutral and unscreened interactions are calculated, the net charge is -0.000357628
notice(2): charge.pppm: RMS error: 0.000406376 Time 00:02:24 | Step 8 / 250 | TPS 0.770384 | ETA 00:05:14 Time 00:02:35 | Step 19 / 250 | TPS 1.03705 | ETA 00:03:42 Time 00:02:45 | Step 30 / 250 | TPS 1.07384 | ETA 00:03:24 Time 00:02:56 | Step 41 / 250 | TPS 1.0386 | ETA 00:03:21 Time 00:03:07 | Step 52 / 250 | TPS 1.03282 | ETA 00:03:11 Time 00:03:17 | Step 63 / 250 | TPS 1.04259 | ETA 00:02:59 Time 00:03:28 | Step 74 / 250 | TPS 1.04234 | ETA 00:02:48 Time 00:03:39 | Step 85 / 250 | TPS 1.00901 | ETA 00:02:43 Time 00:03:49 | Step 96 / 250 | TPS 1.04988 | ETA 00:02:26 Time 00:03:59 | Step 107 / 250 | TPS 1.06528 | ETA 00:02:14 Time 00:04:10 | Step 118 / 250 | TPS 1.03596 | ETA 00:02:07 Time 00:04:21 | Step 129 / 250 | TPS 1.03033 | ETA 00:01:57 Time 00:04:31 | Step 140 / 250 | TPS 1.0352 | ETA 00:01:46 Time 00:04:42 | Step 151 / 250 | TPS 1.03773 | ETA 00:01:35 Time 00:04:52 | Step 162 / 250 | TPS 1.0362 | ETA 00:01:24 Time 00:05:03 | Step 173 / 250 | TPS 1.02768 | ETA 00:01:14 Time 00:05:14 | Step 184 / 250 | TPS 1.0412 | ETA 00:01:03 Time 00:05:24 | Step 195 / 250 | TPS 1.03206 | ETA 00:00:53 Time 00:05:35 | Step 206 / 250 | TPS 1.05436 | ETA 00:00:41 Time 00:05:46 | Step 217 / 250 | TPS 1.02297 | ETA 00:00:32 Time 00:05:56 | Step 228 / 250 | TPS 1.03952 | ETA 00:00:21 Time 00:06:07 | Step 239 / 250 | TPS 1.04558 | ETA 00:00:10 Time 00:06:17 | Step 250 / 250 | TPS 1.0911 | ETA 00:00:00 Average TPS: 1.03 --------- -- Neighborlist stats: 64 normal updates / 0 forced updates / 0 dangerous updates n_neigh_min: 235 / n_neigh_max: 561 / n_neigh_avg: 420.825 shortest rebuild period: 3 -- Cell list stats: Dimension: 3, 3, 3 n_min : 175 / n_max: 292 / n_avg: 222.222 ** run complete **
Now we load the RDF we computed and plot it. Note this is not actually coarse-grained, we're just looking at the mapped rdf which can then be used for coarse-graining.
## Plot RDFs between Mapped CG beads (center of masses) and the reference (C-O in all atom system )
matplotlib.style.use('seaborn-muted')
aa_rdf = model.avg_aa_rdf.result().numpy()
cg_rdf = model.avg_cg_rdf.result().numpy()
r = aa_rdf[1,:]
plt.plot(aa_rdf[1,:], aa_rdf[0,:] / aa_rdf[0,-1],label ='All Atom C-C')
plt.plot(cg_rdf[1,:], cg_rdf[0,:] / cg_rdf[0,-1], label='Mapped (CG)')
plt.xlim(2.5,9)
plt.ylim(0,3)
plt.xlabel(r'r [$\AA$]')
plt.ylabel('$g(r)$')
plt.legend()
plt.show()