making rigid bodies -beginner - please help

[Bojana Zivkovic]

im new to hoomd so this question is most likely stupid but i need help..

im having two types of particles: nanoparticles and ligands.. nanoparticles should stay as just one bead but i need to make a rigid bodies out of ligands..

my central particle should be the same as constituents , since im making the rigid body of just one type of particles(ligands)?? and in my code i have two constituents- ig i wanted for example to make my rigid body to look like aromatic ring or longer alkane chain i just add more const.particles???

this is the error that i get when i run my code: Error intializing ForceComposite

from __future__ import division

import hoomd

from hoomd import md

from hoomd import deprecated

import gsd.hoomd

import random

import numpy as np

import math

# creating a box

        

def hoomdRun():

    sig = [2.0, 0.5]

    #rc = max(sig)*2.5

    delta = [(sig[0]+sig[0])/2 - 1, (sig[1]+sig[1])/2 - 1]

    ijrange = 2

    types = ['1','2']

    mass = [1,0.5]

    initfile = 'nano_ligand_4.gsd'

    log_name = 'output_file'

    finalxml = 'final.xml'

    finaldcd = 'final.dcd'

    movie = 1000     #how often(timesteps) we want to make a movie

    dT = 0.0025   #timestep

    eps = [1,0.5] #need the same number of epsilons as sigmas

    diameter = [sig[0],sig[1]]

    typeid = [0,1]

    rcD = [sig[0]*2**(1/6) + delta[0],((sig[1]+sig[0])/2.0)*2.5,sig[1]*2**(1/6) + delta[1]]   # now our rc = rc + delta for SLJ.. for 12 interactions delta is 0 and for 11 and 22 delta  (sig1+sig1)/2 - 1 or (sig2+sig2)/2 - 1

    

    #creating box with 900atoms

    

    m = 10 #molecules/length

    a = 3 #distance between molecule

    

    box_size = m*a + a

    coord = []

    #initialize x low, y low, z low to build 3d grid.....setting the positions for molecules(nanoparticles and ligands)

    xlo = -a*m/2.0 + a/2.0

    ylo = -a*m/2.0 + a/2.0

    zlo = -a*m/2.0 + a/2.0

    for i in xrange(m): # z

        zlo += a # moving up the z

        for j in xrange(m): # y

            ylo +=a # moving up the y

            

            for k in xrange(m): #x

                coord.append([xlo,ylo,zlo]) #x,y,z coords

                xlo += a

    

            xlo = -a*m/2.0 + a/2.0 #resets xlo

        ylo = -a*m/2.0 + a/2.0 # resets ylo

    #assigning the paramters for each molecule of the system

    num_types = [111,889]  # number of nanoparticles and ligands respectively

    

    diameter_list, types_id_list, mass_list = atom_gsd(diameter,typeid, mass, num_types)

    # create a snapshot

    s = gsd.hoomd.Snapshot()

    s.particles.N = num_types[0]+num_types[1]            # number of rigid bodies

    numP = s.particles.N     #number of particles

    s.particles.types = [types[0],types[1]]

    s.particles.typeid = types_id_list[:]

    s.particles.position = coord[:]

    # s.particles.body = [0,1,2,3]         #body id 4 diff bodies but they can have same types and typeide used for rigid body

    s.particles.mass = mass_list[:]           # mass of nanopart and ligands

    s.particles.charge = numP*[0]

    s.particles.diameter = diameter_list[:]  # diamater is our sigma for diff part

    

    # s.particles.moment_inertia = [numP*[Idiag]]    not need it for this run but i will need it in future

    s.configuration.box = [box_size,box_size,box_size,0,0,0]

    s.configuration.dimensions = 3

    gsd.hoomd.create(name=initfile,snapshot = s)

    #running hoomd

    hoomd.context.initialize("")

    system = hoomd.init.read_gsd(filename=initfile)

    nl = md.nlist.cell()  #creating neighbour list

    #rbuff, checkperiod = nl.tune(warmup=200000, r_min=0.05, r_max=1.0, jumps=20, steps=5000, set_max_check_period=False)

    rigid = hoomd.md.constrain.rigid()

    system.particles.types.add('2')

    rigid.set_param('2',positions=[(0.5,0,0)],types=['2'], diameters=[sig[0],sig[1]])

    rigid.enable()

    rigid.create_bodies()

    center = hoomd.group.rigid.center() ##this is when we have a rigid body but for right now we just have particles

    part = hoomd.group.all() #all the particles and their attributs like position etc

    count = 0  #using count so when for loop goes through 11 12 22 it adds it automatically

    slj = md.pair.slj(r_cut= max(rcD),nlist=nl,d_max=max(sig)) #used max[rc] so we dont have to figure out on our own which one is the max cutoff even tho for this case it is easy to figure it out

    for i in range(0,ijrange):

        for j in range(i,ijrange):

            slj.pair_coeff.set(types[i],types[j],epsilon=math.sqrt(eps[i]*eps[j]),sigma=0.5*(sig[i]+sig[j]),r_cut=rcD[count]) ###setting pair potentials and connect them with the types(11,12,22)

            count += 1

    slj.set_params(mode="shift")

    ## enforces 2d ---- hoomd.md.update.enforce2d()

    hoomd.md.integrate.mode_standard(dt=dT) #setting up the integration for the timestep we choose

    rNUM = random.randint(1,100000)   #random seed generater , important for LV thermostat

    lang = hoomd.md.integrate.langevin(group=part, kT=0.2, seed = rNUM)

    logger = hoomd.analyze.log(filename=log_name, quantities=['temperature','kinetic_energy','potential_energy','pressure','volume','total_energy'],period=movie,header_prefix='#',overwrite=True)

    logger.register_callback('total_energy' , lambda TE : logger.query('kinetic_energy')+logger.query('potential_energy'))

    hoomd.deprecated.dump.xml(group=part,filename =finalxml,vis=True,image=True)#period=100)  vis=True helping with periodic boundary conditions

    hoomd.dump.dcd(filename=finaldcd,period=movie, group = part ,overwrite = True) #owerwrite true -everytime you run the simulation it owerwrites it..False would append to the previous one

    hoomd.run(100000)

    lang.disable()

    nve =hoomd.md.integrate.nve(group=part)

    hoomd.run(100000)

def atom_gsd(diameter,typeid, mass, num_types):

    

    diameter_list = []

    types_id_list = []

    mass_list = []

    

    for i in xrange(len(num_types)):

        for j in xrange(num_types[i]):

            diameter_list.append(diameter[i])

            types_id_list.append(typeid[i])

            mass_list.append(mass[i])

    

    """print 'diamater_list_length ', len(diameter_list)

    print 'types_id_list_lenght ', len(types_id_list)

    print 'types_list_lenght ', len(types_list)

    print 'mass_list ' , len(mass_list)"""

    return diameter_list, types_id_list, mass_list

hoomdRun()

[Joshua Anderson]

The actual error is printed a few lines up from the one you quote: "A rigid body type may not contain constituent particles that are also rigid bodies!".

In hoomd, you must define particles at the center of mass of the rigid bodies with a type that is unique only to those central particles. If you have an actual interaction site at the center of mass, it still needs to be a unique type, but may have identical pair interaction parameters. If you have no interaction sites at the center of mass, then you need to define a type that does not interact with other particles.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan
http://www-personal.umich.edu/~joaander/

--
You received this message because you are subscribed to the Google Groups "hoomd-users" group.
To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users+unsubscribe@googlegroups.com.
To post to this group, send email to hoomd...@googlegroups.com.
Visit this group at https://groups.google.com/group/hoomd-users.
For more options, visit https://groups.google.com/d/optout.

[Nathan Horst]

Seeing this message, I am still a bit confused as to how central particles work exactly. If I define a central particle at the center of mass of a rigid body (let's say at [0,0,0]) and I would like a bond connecting this particle to another rigid body center (at [1,1,1]), does this mean that I would need to define both a central particle and a constituent particle at [0,0,0] and two more at [1,1,1]?? As i understand from the message above, the central particle should not have any explicit pairwise or bonded interactions, and this should be handled by a constituent? If this is the case, do we still have to define pairwise interactions for the central particle type as well, or will HOOMD know that no such interactions exist? Perhaps I am understanding your above message incorrectly.

Thanks,

Nathan Horst

To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users...@googlegroups.com.

[Joshua Anderson]

Reread my statement again, I carefully classified each of the clauses ("If you have an actual interaction site at the center of mass", ...  "If you have no interaction sites at the center of mass").

To answer your specific question, the central particles are real actual particles as far as all of HOOMD is concerned, they may have bonds attached to them and they must have pair coefficients assigned to their types. Those pair coefficients may be zero or non-zero. The constituent particles are also real particles and may have bonds attached to them, and must have pair coefficients assigned. I cannot think of any use-cases where you would need two particles at the same position in space.

I will write a tutorial to demonstrate how to use constrain.rigid. Would these examples be sufficient to clarify the varying use-cases for everyone?

1. Rigid rods: 9 particles long. Demonstrate that the central particle may have identical interactions as constituents, and how hoomd can create constituent particles.

2. A short polymer tethered to the corner of a rigid cube. Demonstrate that the central particle may be non-interacting and that bonds may attach to constituent particles, and how a user can specify the constituent particles.

------
Joshua A. Anderson, Ph.D.
Research Area Specialist, Chemical Engineering, University of Michigan
http://www-personal.umich.edu/~joaander/

To unsubscribe from this group and stop receiving emails from it, send an email to hoomd-users+unsubscribe@googlegroups.com.

[Nathan Horst]

That tutorial seems like it would be a great supplement to the current documentation on constrain.rigid. 

Thanks for the clarification,

Nathan

[Xinyan Yang]

Are these tutorials available somewhere? I could not find them in the HOOMD documentation.

Best,

Xinyan

[Joshua Anderson]

Xinyan,

You can find the current rigid body tutorial here, linked to from the documentation: https://hoomd-blue.readthedocs.io/en/v4.3.0/tutorial/06-Modelling-Rigid-Bodies/00-index.html

If you want to find the rigid body tutorial as it existed in this message you replied to from 6 years ago, then you can examine the git repository log and check out the appropriate commit.

------

Joshua A. Anderson, Ph.D.

Research Area Specialist, Chemical Engineering, University of Michigan

To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/56a06c92-903d-4260-86d0-060ba86a2379n%40googlegroups.com.

[Xinyan Yang]

Great! Thank you!

Best,

To view this discussion on the web visit https://groups.google.com/d/msgid/hoomd-users/41C2728D-46AF-4AB1-8366-75052908D603%40umich.edu.