On Thu, Jun 29, 2017 at 7:08 PM, Thomas Evangelidis <tevang3.gmail.com>
wrote:
> Greetings,
>
> I am using the Parmed API to prepare a system for simulation with OpenMM. I
> want to add a restraint 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 (or a
> virtual site maybe?) to the system with mass=0 and fixed position. However,
>
You can't use a mass of 0 with restraints -- you'll get an infinite
acceleration (it'll work with constraints). Also, the smaller the mass you
use, the shorter the time step will need to be to deal with the large
resulting accelerations and velocities.
> I still couldn't find how to set the dummy atom's coordinates to x0,y0,z0.
> This is how my code looks like:
>
>
>
> from parmed.amber import AmberParm
>
import parmed as pmd
# Needed to add an atom
import copy # To create a dummy atom
> # Load the Amber files
> parm = AmberParm(args.PRMTOP, args.COORD)
This is where you should add the atom. So you can do something like this:
parm.add_atom(copy.copy(parm.atoms[-1]), 'DUM', parm.residues[-1].number +
1)
# Now set various attributes of the atom to make it look different
dummy = parm.atoms[-1]
dummy.mass = 12.01
dummy.charge = 0
dummy.name = 'DU'
dummy.xx, dummy.xy, dummy.xz = x0, y0, z0
parm.remake_parm() # Needed to update the data that actually gets written
to the prmtop file
# Give the dummy atom 0 Lennard-Jones
pmd.actions.addLJType(parm, '.DU', radius=1.0, epsilon=0.0).execute()
> # 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 also the function
> parm.add_atom which seemingly does the same thing. Which is the right way
> to do it?
>
What you did above was to create an OpenMM System object. This is used to
run an MD simulation on GPUs using OpenMM (instead of Amber), but using the
Amber force field. If your intention is to use OpenMM for the MD
simulation, then this is the next step (and consult the OpenMM docs for how
to run a simulation).
If all you want to do is add the dummy atom and generate a prmtop and
inpcrd you can use with Amber, delete the system creation part. Either
way, "add_atom" is what you want to use.
If you want to save a prmtop and inpcrd, just do this:
parm.save('prmtop_with_dummy.parm7') # .prmtop suffix also works
parm.save('inpcrd_with_dummy_atom.rst7') # Writes the inpcrd file
HTH,
Jason
--
Jason M. Swails
_______________________________________________
AMBER mailing list
AMBER.ambermd.org
http://lists.ambermd.org/mailman/listinfo/amber
Received on Thu Jun 29 2017 - 19:00:02 PDT
