Page 1 of 1

adding dummy atom

Posted: Thu Jun 29, 2017 3:51 pm
by tevang
Greetings,

I want to add a CustomCentroidBondForce between the COM of the ligand and it's initial COM coordinates x0,y0,z0. For this purpose I thought to add a dummy atom to the system with mass=0. However, I still can't find how to set it's coordinates to x0,y0,z0. This is how my code looks like:

Code: Select all

# Load the Amber files
parm = AmberParm(args.PRMTOP, args.COORD)

# Create the OpenMM system
# Load the Amber topology file
system = parm.createSystem(nonbondedMethod=CutoffNonPeriodic,
                nonbondedCutoff=args.NONBONDED_CUTOFF*u.angstrom,
                constraints=app.HBonds, rigidWater=True,
                implicitSolvent=igb,
                # implicitSolventKappa=0.107317448721*(1.0/u.angstrom),   # the value for 0.2 M salt
                implicitSolventSaltConc=args.SALTCON*u.moles/u.liter,
                soluteDielectric=1.0,
                solventDielectric=78.5,
                removeCMMotion=True,
                ewaldErrorTolerance=5e-05,
                flexibleConstraints=False,
                verbose=args.VERBOSE,
                useSASA=False
    
COM_index = system.addParticle(0.0)    
And then..? How do I set the coordinates? There is another function named parm.add_atom, but it is relevant to AmberParm. Which is the right way to do it?

Re: adding dummy atom

Posted: Thu Jun 29, 2017 4:16 pm
by peastman
A much easier way is to use per-bond parameters for the coordinates:

Code: Select all

force = CustomCentroidBondForce(1, "k*((x1-x0)^2+(y1-y0)^2+(z1-z0)^2)")
force.addPerBondParameter("x0")
force.addPerBondParameter("y0")
force.addPerBondParameter("z0")
force.addGroup(particles)
force.addBond([0], [x0, y0, z0])
Peter

Re: adding dummy atom

Posted: Thu Jun 29, 2017 5:30 pm
by tevang
You are right. Damn it! This is exactly what I did at the beginning but I had

Code: Select all

force.addPerBondParameter('x0')
force.addPerBondParameter('y0')
force.addPerBondParameter('z0')
force.addPerBondParameter("k")
force.addGroup(ligatoms)
force.addBond([0], [k, x0, y0, z0])
and thus the ligand was drifting.