Neibourlist issue

[Shailesh-sathe]

import hoomd
import gsd.hoomd
import hoomd.snapshot
import hoomd.state
import numpy as np
import os

N =500
sigma = 0.5
epsilon = 1
KT = 1.0
L=17.5
runs=5e4
buffers=[0.2,0.3,0.4,0.6,0.8,1,2,3,5]
for Buffer in buffers:
os.makedirs('checktraj',exist_ok=True)
output=f'checktraj/check4_{Buffer}.gsd'
if os.path.exists(output):
os.remove(output)

class writingframes(hoomd.custom.Action):
    def __init__(self, filename):
        self.filename = filename
        self.gsd_file = gsd.hoomd.open(name=filename, mode='w')
        self.frames_written = 0

    def act(self, timestep):
        snapshot = self._state.get_snapshot()
        if snapshot.communicator.rank == 0:
            frame = gsd.hoomd.Frame()
            frame.configuration.step = timestep
            frame.configuration.box = snapshot.configuration.box
            frame.particles.N = snapshot.particles.N
            frame.particles.types = snapshot.particles.types
            frame.particles.typeid = snapshot.particles.typeid
            frame.particles.position = snapshot.particles.position
            frame.bonds.N = snapshot.bonds.N
            frame.bonds.types = snapshot.bonds.types
            frame.bonds.typeid = snapshot.bonds.typeid
            frame.bonds.group = snapshot.bonds.group
            frame.particles.velocity=snapshot.particles.velocity
            # Validate frame before writing
            assert frame.particles.position.shape[0] == frame.particles.N, \
                f"Position count {frame.particles.position.shape[0]} does not match particle count {frame.particles.N}"
            assert frame.bonds.group.shape[0] == frame.bonds.N, \
                f"Bond group count {frame.bonds.group.shape[0]} does not match bond count {frame.bonds.N}"

            # Write the frame to the GSD file
            self.gsd_file.append(frame)
            self.frames_written += 1
            print(f"Frame {self.frames_written} written at timestep {timestep}")

class BondTuner(hoomd.custom.Action):
    def __init__(self,bond_force,initK,finalK,rate):
        self.bond_force=bond_force
        self.k=initK
        self.finalK=finalK
        self.rate=rate
    def act(self, timestep):
        snapshot = self._state.get_snapshot()
        if self.k < self.finalK:
            self.k += self.rate
            self.bond_force.params['Crosslink'] = dict(k=self.k, r0=1.0)
            if(snapshot.communicator.rank==0):
                print(f"At timestep {timestep}, updated bond k to {self.k}")

# class softtuner(hoomd.custom.Action):
#     def __init__(self,pair_force,initA,finalA,rate):
#         self.pair_force=pair_force
#         self.k=initA
#         self.finalK=finalA
#         self.rate=rate
#     def act(self, timestep):
#         if self.k < self.finalK:
#             self.k += self.rate
#             self.pair_force.params[('A', 'A')] = dict(A=self.k, rho=2, C=0)
#             self.pair_force.params[('A', 'B')] = dict(A=self.k, rho=2, C=0)
#             self.pair_force.params[('A', 'C')] = dict(A=self.k, rho=2 , C=0)
#             self.pair_force.params[('B', 'B')] = dict(A=self.k, rho=2, C=0)
#             self.pair_force.params[('B', 'C')] = dict(A=self.k, rho=2, C=0)
#             self.pair_force.params[('C', 'C')] = dict(A=self.k, rho=2, C=0)
#             print(f"At timestep {timestep}, updated bond k to {self.k}")

class HarmonicCylinder(hoomd.md.force.Custom):
    def __init__(self, radius, k, axis=(0, 0, 1)):
        # Initialize the CustomForceCompute with the provided trigger
        super().__init__()
        
        self.radius = radius
        self.k = k
        self.axis = np.array(axis) / np.linalg.norm(axis)

    def set_forces(self, timestep):
        # print(timestep)
        i=0
        with self.cpu_local_force_arrays as arrays, self._state.cpu_local_snapshot as snap:
            forces = np.zeros_like(snap.particles.position)
            positions = snap.particles.position
            nparticles=len(positions)
            for  i in range(nparticles):
                force_array = arrays.force
                radial_pos=np.zeros(3)
                radial_pos[0:2] = positions[i][0:2]  # Project onto plane perpendicular to axis
                radial_dist = np.linalg.norm(radial_pos)
                # print(radial_dist,positions[i])
                # print(pos)
                if radial_dist > self.radius:
                    # print(i) 
                    displacement = radial_dist - self.radius
                    radial_dir = radial_pos / radial_dist
                    force_magnitude =   - self.k * displacement
                    forces[i] = force_magnitude * radial_dir
                    force_array[i] = forces[i]

class StepWriter(hoomd.custom.Action):

    def act(self, timestep):
        snap=self._state.get_snapshot()
        if snap.communicator.rank==0:
            print(timestep)

def main(buffer):
    N =500
    sigma = 0.7
    epsilon = 1
    KT = 1.0
    L=17.5
    A=5
    runs=5e4
    DD = hoomd.device.auto_select()
    randmseed=np.random.randint(1,1000)
    sim = hoomd.Simulation(device=DD, seed=3)

    # setting harmonic potential
    harmonic = hoomd.md.bond.Harmonic()
    harmonic.params['A-A'] = dict(k=100, r0=1.0)
    harmonic.params['Crosslink'] = dict(k=100, r0=1.0)

    # setting LJ potential
    
    # nlist = hoomd.md.nlist.Stencil(cell_width=10,buffer=0.8, rebuild_check_delay=1, check_dist=False,default_r_cut=2)
    # nlist=hoomd.md.nlist.Tree(buffer=1,rebuild_check_delay=1,check_dist=True,default_r_cut=5)
    nlist=hoomd.md.nlist.Cell(buffer=buffer,rebuild_check_delay=1,check_dist=False,default_r_cut=2)
    
    lj = hoomd.md.pair.LJ(nlist=nlist)
    lj.params[('A', 'A')] = dict(epsilon=1, sigma=sigma)
    lj.params[('B', 'B')] = dict(epsilon=1, sigma=sigma)  # dict(epsilon=0,sigma=0)
    lj.params[('C', 'C')] = dict(epsilon=0, sigma=sigma)
    lj.params[('A', 'B')] = dict(epsilon=1, sigma=sigma)  # dict(epsilon=0,sigma=0)
    lj.params[('A', 'C')] = dict(epsilon=0, sigma=sigma)
    lj.params[('B', 'C')] = dict(epsilon=0, sigma=sigma)  # dict(epsilon=0,sigma=0)
    lj.r_cut[('A', 'A')] = (2.0 ** (1 / 6)) * sigma
    lj.r_cut[('B', 'B')] = (2.0 ** (1 / 6)) * sigma
    lj.r_cut[('C', 'C')] = (2.0 ** (1 / 6)) * sigma
    lj.r_cut[('A', 'B')] = (2.0 ** (1 / 6)) * sigma
    lj.r_cut[('A', 'C')] = (2.0 ** (1 / 6)) * sigma
    lj.r_cut[('B', 'C')] = (2.0 ** (1 / 6)) * sigma


    #Creating soft potential for replicating monomers
    soft=hoomd.md.pair.Buckingham(nlist,default_r_cut=0.8)
    soft.params[('A', 'A')] = dict(A=A, rho=0.6, C=0)
    soft.params[('A', 'B')] = dict(A=A, rho=0.6, C=0)
    soft.params[('A', 'C')] = dict(A=A, rho=0.6 , C=0)
    soft.params[('B', 'B')] = dict(A=A, rho=0.6, C=0)
    soft.params[('B', 'C')] = dict(A=A, rho=0.6, C=0)
    soft.params[('C', 'C')] = dict(A=A, rho=0.6, C=0)

    # creating hard box
    z=L/2
    halfsigma=0.4
    top = hoomd.wall.Plane(origin=(0, z, 0), normal=(0, -1, 0))
    bottom = hoomd.wall.Plane(origin=(0, -z, 0), normal=(0, 1, 0))
    left = hoomd.wall.Plane(origin=(-z, 0, 0), normal=(1, 0, 0))
    right = hoomd.wall.Plane(origin=(z, 0, 0), normal=(-1, 0, 0))
    front = hoomd.wall.Plane(origin=(0, 0, z), normal=(0, 0, -1))
    back = hoomd.wall.Plane(origin=(0, 0, -z), normal=(0, 0, 1))
    ljbox = hoomd.md.external.wall.LJ([top, bottom, left, right, back, front])
    ljbox.params['A'] = dict(epsilon=1, sigma=halfsigma, r_cut=(2.0 ** (1 / 6)) * halfsigma)
    ljbox.params['B'] = dict(epsilon=1, sigma=halfsigma, r_cut=(2.0 ** (1 / 6)) * halfsigma)
    ljbox.params['C'] = dict(epsilon=1, sigma=halfsigma, r_cut=(2.0 ** (1 / 6)) * halfsigma)

    #Creating cylindrical confinement
    cylinder=hoomd.wall.Cylinder(radius=3.50,axis=(0,0,1))
    ljcylinder=hoomd.md.external.wall.LJ([cylinder])
    ljcylinder.params['A'] = dict(epsilon=1, sigma=halfsigma, r_cut=(2.0 ** (1 / 6)) * halfsigma)
    ljcylinder.params['B'] = dict(epsilon=1, sigma=halfsigma, r_cut=(2.0 ** (1 / 6)) * halfsigma)
    ljcylinder.params['C'] = dict(epsilon=1, sigma=halfsigma, r_cut=(2.0 ** (1 / 6)) * halfsigma)

    # setting langevin

    lang = hoomd.md.methods.Langevin(filter=hoomd.filter.All(), kT=KT)
    # lang = hoomd.md.methods.ConstantVolume(filter=hoomd.filter.All(),thermostat=hoomd.md.methods.thermostats.Bussi(kT=1))

    integrator = hoomd.md.Integrator(dt=0.01, methods=[lang])
    cylinder=HarmonicCylinder(radius=1.5,k=100)
    # integrator.forces.append(lj)
    integrator.forces.append(ljcylinder)
    integrator.forces.append(ljbox)
    # integrator.forces.append(harmonic)
    integrator.forces.append(soft)
    # integrator.forces.append(cylinder)
    sim.operations.integrator = integrator

    sim.create_state_from_gsd('SPequiled_3.50_10_0.6.gsd')
    # sim.state.thermalize_particle_momenta(filter=hoomd.filter.All(), kT=1)
    # sim.run(0)
    sim.operations.integrator = integrator
    # writer2 = hoomd.write.GSD(filename=f'{N}polymertrajectoryidealring{i}.gsd', trigger=hoomd.trigger.Periodic(5000), mode='wb')
    writer1=hoomd.write.CustomWriter(trigger=hoomd.trigger.Periodic(1000),action=writingframes(output))
    sim.operations.writers.append(writer1)
    # sim.operations.writers.append(writer2)
    # bondtuner=BondTuner(harmonic,1,100,0.01)
    # custombond=hoomd.update.CustomUpdater(trigger=hoomd.trigger.Periodic(10),action=bondtuner)
    # sim.operations.updaters.append(custombond)
    # step=hoomd.write.CustomWriter(trigger=hoomd.trigger.Periodic(1000),action=StepWriter())
    # sim.operations.writers.append(step)
    if DD.communicator.rank==0:
        print('Run Start')
    sim.operations.computes.append(nlist)
    sim.run(1e5)
    if DD.communicator.rank==0:
        print(f"Initialized nlist with buffer={nlist.cpu_local_nlist_arrays}, check_dist={nlist.check_dist}, default_r_cut")
    if DD.communicator.rank==0:
        print("Run end")
    # sim.operations.integrator.forces.append(lj)
    # sim.operations.integrator.forces.remove(cylinder)
    # sim.operations.integrator.forces.append(ljcylinder)
    # sim.operations.integrator.forces.append(soft)
    # customsoft=hoomd.update.CustomUpdater(trigger=hoomd.trigger.Periodic(100),action=softtuner(soft,0,100,0.5))
    # sim.operations.updaters.append(customsoft)
    # print('Run Start')
    # sim.run(5e5)
    # print("Run end")   
    # # sim.operations.integrator.forces.remove(soft)
    # # sim.operations.integrator.forces.append(lj)
    # sim.run(5e4)
    # writer2.flush()
    # # Open each trajectory file in read mode
    # with gsd.hoomd.open(f'{N}polymertrajectoryidealring{i}.gsd', mode='r') as traj:
    #     for frame in traj:
    #         # Append each frame to the output file
    #         output.append(frame)
    # # Delete the trajectory file after merging
    # os.remove(f'{N}polymertrajectoryidealring{i}.gsd')
    fsnapshot=sim.state.get_snapshot()
    if fsnapshot.communicator.rank == 0:
        # frame = gsd.hoomd.Frame()
        # frame.configuration.box = [40,40,40,0,0,0]
        # frame.particles.N = fsnapshot.particles.N
        # frame.particles.types = fsnapshot.particles.types
        # frame.particles.typeid = fsnapshot.particles.typeid
        # frame.particles.position = fsnapshot.particles.position
        # frame.bonds.N = fsnapshot.bonds.N
        # frame.bonds.types = ['A-A', 'Crosslink']
        # frame.bonds.typeid = [0]*506
        # frame.bonds.group = fsnapshot.bonds.group
        
        # equiledfile = gsd.hoomd.open(f'mymodel2equiled_3.50_{A}_0.6.gsd', mode='w')
        # equiledfile.append(frame)
        print("Snapshot successfully written to the GSD file.")

main(Buffer)

This is my code as I change the neighbourlist parameter my trajectory changes. even though I have set check_dist=False still by changing buffer it changes trajectory and speed of simulation

[joaander]

Yes, I expect changing the value of buffer to change the speed of the simulation. I also expect it to change the specific trajectory followed. When you change buffer, you will change the order of the particles in the neighbor list. Floating point addition is not associative, therefore this leads to different trajectories. Parallel execution on the GPU will also produce different trajectories - even when all parameters are kept identical between two simulation runs.

[Shailesh-sathe]

Thanks for your time sir, but I still don't get the point if trajectories changes then which trajectory is true trajectory?

Whats the point of check_dist=False? it should rebuild the neighbourlist irrespective of buffer. but it's not rebuilding neighbourlist.
Also if trajectories should be same if I am not using GPU and mpi if neighbourlist check_dist=False and I am changing buffer.
Does the range of calculating nlist depend on buffer or only depends on default_r_cut?

[joaander]

Thanks for your time sir, but I still don't get the point if trajectories changes then which trajectory is true trajectory?

MD is an approximate solution to a system of differential equations. In most MD simulations, an analytical solution is not possible so the "true" trajectory is an unknowable quantity. You can learn more about this topic in MD textbooks such as Allen and Tildesley or Frenkel and Smit

Whats the point of check_dist=False?

Some MD codes default to building neighbor lists on a fixed schedule, which is aggressive and can lead to forces being incorrectly computed as 0. check_dist=False exists in HOOMD so that users can fairly compare simulation performance between HOOMD and such codes that do not correctly check the distance. You should always use check_dist=True (the default) in HOOMD to ensure that all pair forces are correctly computed.

it should rebuild the neighbourlist irrespective of buffer. but it's not rebuilding neighbourlist.

It is rebuilding the neighbor list:

hoomd-blue/hoomd/md/NeighborList.cc

Lines 1376 to 1387 in 24ec50e

	if (m_r_buff < 1e-6
	|| (!m_dist_check
	&& (m_rebuild_check_delay == 0
	|| (m_rebuild_check_delay > 1
	&& timestep == (m_last_updated_tstep + m_rebuild_check_delay)))))
	{
	result = true;
	}
	else
	{
	result = distanceCheck(timestep);
	}

Also if trajectories should be same if I am not using GPU and mpi if neighbourlist check_dist=False and I am changing buffer.

In general, you should assume that every time you run a simulation - even with identical parameters - that the resulting trajectory will be different. There may be some rare cases where you can obtain an identical trajectory, but I cannot possibly enumerate all of the special conditions necessary to achieve that.

Does the range of calculating nlist depend on buffer or only depends on default_r_cut?

You should leave the neighbor list's default_r_cut set to 0 when using it in conjunction with a pair potential. The neighbor list will automatically use the maximum r_cut value from all associated pair potentials.

[Shailesh-sathe]

Thank you sir.
What should be buffer size if I am simulating polymer in cylindrical confinement?

[joaander]

Follow this tutorial to find out: https://hoomd-blue.readthedocs.io/en/v4.8.2/howto/choose-the-neighbor-list-buffer-distance.html

[Shailesh-sathe]

Thank you