Continuing with the question from UAMMD

[biubiuubiuu]

In the issue mentioned in the UAMMD repository (https://github.com/RaulPPelaez/UAMMD/issues/38), I followed your advice to modify the JSON file and attempted to remove the extra forces, but the system still "exploded." I then tried reducing the particle count to 12 while keeping the box size the same, and after adding the parameter file in the force field, the simulation ran successfully. Could it be that the particles near the boundaries are experiencing unreasonable forces, leading to the system errors? So, should I regenerate my system, or do you have any better suggestions? Thank you for your patience.

[PabloIbannez]

I have added the bugfix to UAMMD-structured. You need to recompile UAMMD-structured (from scratch) and it should work. I have also added a PSE example to the Examples folder.

Feel free to share the code you use to generate the .json input file here and I give you a hand to simulate the system you are interested in.

[biubiuubiuu]

I have added the bugfix to UAMMD-structured. You need to recompile UAMMD-structured (from scratch) and it should work. I have also added a PSE example to the Examples folder.

Feel free to share the code you use to generate the .json input file here and I give you a hand to simulate the system you are interested in.

Thank you for your timely update. In my system, I process a protein into 12 beads, with bonds and angles defined between the beads. The system contains multiple proteins. However, I understand that PSE treats the whole system as a colloidal particle. Does this mean that the interactions between the beads within a protein will be ignored, and only hydrodynamic interactions will be considered?

[PabloIbannez]

The interaction is always between the particles. Here you have an example of a system I think could be similar to yours:

import pyUAMMD

import numpy as np

def createProtein(n, nPartPerProtein, L, radius):

    firstParticlePos = np.random.uniform(-0.5*L, 0.5*L, 3)

    firstParticleVector = np.random.uniform(-1, 1, 3)
    firstParticleVector /= np.linalg.norm(firstParticleVector)

    # You can write here the model of your protein
    # For example, a simple linear protein (like a polymer, using the Worm Like Chain model)

    # We create a protein of len nPartPerProtein
    # The position of the first particle is firstParticlePos
    # and the direction of the protein is firstParticleVector.
    # The distance between particles is 2*radius

    # The ids of protein n start at n*nPartPerProtein
    idOffset = n*nPartPerProtein
    ids = np.arange(nPartPerProtein)
    ids = ids + idOffset

    # Convert ids to list of regular integers
    ids = ids.tolist()
    ids = [int(i) for i in ids]

    particles = []
    for i in range(nPartPerProtein):
        particles.append(firstParticlePos + i*2*radius*firstParticleVector)

    bondList = []
    for i in range(nPartPerProtein-1):
        bondList.append([idOffset+i, idOffset+i+1])

    angleList = []
    for i in range(nPartPerProtein-2):
        angleList.append([idOffset+i, idOffset+i+1, idOffset+i+2])

    return ids, particles, bondList, angleList

nProteins         = 50
nPartPerProtein   = 12
L                 = 100
radius            = 1.0

Kb = 100.0 # Harmonic bond constant
Ka = 20.0  # Kratky-Porod angular bond constant

numberOfSteps     = 100000
numberStepsOutput = 1000
timeStep          = 0.01

temperature       = 1.0
viscosity         = 1.0
psi               = 0.6
tolerance         = 1e-3

box = [L,L,L]

simulation = pyUAMMD.simulation()

simulation["system"] = {}
simulation["system"]["info"] = {}
simulation["system"]["info"]["type"] = ["Simulation","Information"]
simulation["system"]["info"]["parameters"] = {}
simulation["system"]["info"]["parameters"]["name"] = "ExamplePSE"

simulation["global"] = {}
simulation["global"]["types"] = {}
simulation["global"]["types"]["type"]   = ["Types","Basic"]
simulation["global"]["types"]["labels"] = ["name", "radius", "mass", "charge"]
simulation["global"]["types"]["data"]  = [["A", radius, 0, 0]]

simulation["global"]["ensemble"] = {}
simulation["global"]["ensemble"]["type"]   = ["Ensemble","NVT"]
simulation["global"]["ensemble"]["labels"] = ["box", "temperature"]
simulation["global"]["ensemble"]["data"]   = [[box, temperature]]

simulation["integrator"] = {}
simulation["integrator"]["PSE"] = {}
simulation["integrator"]["PSE"]["type"] = ["BDHITriplePeriodic", "PositivelySplitEwald"]
simulation["integrator"]["PSE"]["parameters"] = {}
simulation["integrator"]["PSE"]["parameters"]["timeStep"]           = timeStep
simulation["integrator"]["PSE"]["parameters"]["viscosity"]          = viscosity
simulation["integrator"]["PSE"]["parameters"]["hydrodynamicRadius"] = radius
simulation["integrator"]["PSE"]["parameters"]["tolerance"]          = tolerance
simulation["integrator"]["PSE"]["parameters"]["psi"]                = psi

simulation["integrator"]["schedule"] = {}
simulation["integrator"]["schedule"]["type"] = ["Schedule", "Integrator"]
simulation["integrator"]["schedule"]["labels"] = ["order", "integrator","steps"]
simulation["integrator"]["schedule"]["data"] = [
    [1, "PSE", numberOfSteps],
]

simulation["state"] = {}
simulation["state"]["labels"] = ["id", "position"]
simulation["state"]["data"] = []

# Positions will be specified after, when proteins are created

simulation["topology"] = {}

simulation["topology"]["structure"] = {}
simulation["topology"]["structure"]["labels"] = ["id", "type"]
simulation["topology"]["structure"]["data"] = []

bondList  = []
angleList = []

for i in range(nProteins):
    ids, particles, bList, aList = createProtein(i, nPartPerProtein, L, radius)

    for index, pid in enumerate(ids):
        simulation["state"]["data"] += [[pid, particles[index].tolist()]]
        simulation["topology"]["structure"]["data"] += [[pid, "A"]]

    bondList  += bList
    angleList += aList

simulation["topology"]["forceField"] = {}
simulation["topology"]["forceField"]["bonds"] = {}
simulation["topology"]["forceField"]["bonds"]["type"] = ["Bond2","Harmonic"]
simulation["topology"]["forceField"]["bonds"]["parameters"] = {}
simulation["topology"]["forceField"]["bonds"]["labels"] = ["id_i", "id_j", "K", "r0"]
simulation["topology"]["forceField"]["bonds"]["data"] = []

for i, bond in enumerate(bondList):
    simulation["topology"]["forceField"]["bonds"]["data"] += [[bond[0], bond[1], Kb, 2*radius]]

simulation["topology"]["forceField"]["angles"] = {}
simulation["topology"]["forceField"]["angles"]["type"] = ["Bond3","KratkyPorod"]
simulation["topology"]["forceField"]["angles"]["parameters"] = {}
simulation["topology"]["forceField"]["angles"]["labels"] = ["id_i", "id_j", "id_k", "K"]
simulation["topology"]["forceField"]["angles"]["data"] = []

for i, angle in enumerate(angleList):
    simulation["topology"]["forceField"]["angles"]["data"] += [[angle[0], angle[1], angle[2], Ka]]

simulation["simulationStep"] = {}
simulation["simulationStep"]["info"] = {}
simulation["simulationStep"]["info"]["type"] = ["UtilsStep", "InfoStep"]
simulation["simulationStep"]["info"]["parameters"] = {"intervalStep":numberStepsOutput}

simulation["simulationStep"]["pos"] = {}
simulation["simulationStep"]["pos"]["type"] = ["WriteStep", "WriteStep"]
simulation["simulationStep"]["pos"]["parameters"] = {"intervalStep":numberStepsOutput,
                                                     "outputFilePath": "output",
                                                     "outputFormat":"xyz"}

simulation.write(f"simulation.json")

This script will generate an input to simulate a system made up of a number of "proteins", where each protein is made up of 12 particles. These particles interact with each other through distance and angle dependent bonds. So the total interaction between the particles would be the hydrodynamic interaction given by PSE plus the internal interactions.