Build a new molecule from scratch by joining the atoms in script

[sajid iqbal]

Hi Sire Community,

I am wondering if Sire facilitates to construct structures of organic molecules by writing script, for example benzene. I know, we can optimize the structures using the Amber .inp and .prm file or PDB files on Sire. Any example in tutorial would be great help.

Regards,

Thanks

[Christopher Woods]

  Dear Sajid,

Thanks for getting in touch and for your question. Yes, it is possible to build a molecule in Sire using a script, although in many cases it is much easier to use the graphical tools that come with Amber. I have attached a commented script that builds a pair of benzene molecules using Sire, and then calculates their intermolecular coulomb and LJ energy. I apologise that the molecule creation API is quite messy and not well documented. It was one of the earliest things written in Sire, but it turned out that we didn't really use it for anything, and building molecules using dedicated tools was easier.

  Best wishes,

  Christopher

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

[Christopher Woods]

Sending again with the script inline, so that it is easier to find using a google search ;-)

=================================

from Sire.Mol import *

from Sire.IO import *

from Sire.Maths import *

from Sire.MM import *

from Sire.Units import *

# create the molecule

benzene = Molecule("benzene")

# create a residue for all of the atoms

res = benzene.edit().add( ResName("BEN") )

res.renumber( ResNum(1) )

# Add a CutGroup - all atoms

# have to be children of a CutGroup (a single entity

# that is used for cutoff or periodic boundary calculations,

# and is essentially the same as a residue) - this requirement

# for cutgroups is very confusing and something I want to remove!

res.molecule().add( CGName("1") )

# Now add the atoms to the residue, ensuring they are 

# a child of the cutgroup

res.add( AtomName("H1") ).reparent( CGName("1") )

res.add( AtomName("H2") ).reparent( CGName("1") )

res.add( AtomName("H3") ).reparent( CGName("1") )

res.add( AtomName("H4") ).reparent( CGName("1") )

res.add( AtomName("H5") ).reparent( CGName("1") )

res.add( AtomName("H6") ).reparent( CGName("1") )

res.add( AtomName("C1") ).reparent( CGName("1") )

res.add( AtomName("C2") ).reparent( CGName("1") )

res.add( AtomName("C3") ).reparent( CGName("1") )

res.add( AtomName("C4") ).reparent( CGName("1") )

res.add( AtomName("C5") ).reparent( CGName("1") )

res.add( AtomName("C6") ).reparent( CGName("1") )

# now we have added all of the atoms, commit this

# so that we can add properties to the molecule

# (properties are the connectivity of the atoms,

#  coordinates of the atoms, charges, LJ parameters etc.)

benzene = res.molecule().commit()

# here are the charges and LJ parameters for benzene

# (please look these up as I have just guessed them!)

c_chg = -0.0622 * mod_electron

h_chg =  0.0622 * mod_electron

c_lj = LJParameter( 3.5 * angstrom,  -0.015 * kcal_per_mol )

h_lj = LJParameter( 1.0 * angstrom, -0.005 * kcal_per_mol )

# let's set the coordinates, element, charges and LJ parameters

benzene = benzene.edit() \

                 .atom( AtomName("C1") ) \

                 .setProperty("coordinates", Vector(0.000,  1.396,  0.000)) \

                 .setProperty("element", Element("C")) \

                 .setProperty("charge", c_chg).setProperty("LJ", c_lj) \

                 .molecule().atom( AtomName("C2") ) \

                 .setProperty("coordinates", Vector(1.209,  0.698,  0.000)) \

                 .setProperty("element", Element("C")) \

                 .setProperty("charge", c_chg).setProperty("LJ", c_lj) \

                 .molecule().atom( AtomName("C3") ) \

                 .setProperty("coordinates", Vector(1.209, -0.698,  0.000)) \

                 .setProperty("element", Element("C")) \

                 .setProperty("charge", c_chg).setProperty("LJ", c_lj) \

                 .molecule().atom( AtomName("C4") ) \

                 .setProperty("coordinates", Vector(0.000, -1.396,  0.000)) \

                 .setProperty("element", Element("C")) \

                 .setProperty("charge", c_chg).setProperty("LJ", c_lj) \

                 .molecule().atom( AtomName("C5") ) \

                 .setProperty("coordinates", Vector(-1.209, -0.698,  0.000)) \

                 .setProperty("element", Element("C")) \

                 .setProperty("charge", c_chg).setProperty("LJ", c_lj) \

                 .molecule().atom( AtomName("C6") ) \

                 .setProperty("coordinates", Vector(-1.209,  0.698,  0.000)) \

                 .setProperty("element", Element("C")) \

                 .setProperty("charge", c_chg).setProperty("LJ", c_lj) \

                 .molecule().atom( AtomName("H1") ) \

                 .setProperty("coordinates", Vector( 0.000,  2.479,  0.000)) \

                 .setProperty("element", Element("H")) \

                 .setProperty("charge", h_chg).setProperty("LJ", h_lj) \

                 .molecule().atom( AtomName("H2") ) \

                 .setProperty("coordinates", Vector( 2.147,  1.240,  0.000)) \

                 .setProperty("element", Element("H")) \

                 .setProperty("charge", h_chg).setProperty("LJ", h_lj) \

                 .molecule().atom( AtomName("H3") ) \

                 .setProperty("coordinates", Vector( 2.147, -1.240,  0.000)) \

                 .setProperty("element", Element("H")) \

                 .setProperty("charge", h_chg).setProperty("LJ", h_lj) \

                 .molecule().atom( AtomName("H4") ) \

                 .setProperty("coordinates", Vector( 0.000, -2.479,  0.000)) \

                 .setProperty("element", Element("H")) \

                 .setProperty("charge", h_chg).setProperty("LJ", h_lj) \

                 .molecule().atom( AtomName("H5") ) \

                 .setProperty("coordinates", Vector(-2.147, -1.240,  0.000)) \

                 .setProperty("element", Element("H")) \

                 .setProperty("charge", h_chg).setProperty("LJ", h_lj) \

                 .molecule().atom( AtomName("H6") ) \

                 .setProperty("coordinates", Vector(-2.147,  1.240,  0.000)) \

                 .setProperty("element", Element("H")) \

                 .setProperty("charge", h_chg).setProperty("LJ", h_lj) \

                 .molecule().commit()

# now set the connectivity of the molecule - this will automatically

# work out the connectivity based on the elements and interatomic distances

connectivity = Connectivity(benzene)

# print it to check it is correct. If not, look at the

# functions of this object to add or remove bonds

print(connectivity)

# set the connectivity property of the molecule

benzene = benzene.edit().setProperty("connectivity", connectivity).commit()

# print out a PDB of the molecule

PDB().write(benzene, "benzene.pdb")

# duplicate the molecule so that we can calculate inter-benzene energies

benzene2 = benzene.edit().renumber().commit()

# translate this second molecule by 5 A along the Z axis

benzene2 = benzene2.move().translate( Vector(0,0,5) ).commit()

# create a forcefield to calculate the coulomb and LJ energy

cljff = InterFF()

cljff.add(benzene)

cljff.add(benzene2)

print("\nCalculating the interaction energy...")

print("Coulomb = %s" % cljff.energy( cljff.components().coulomb() ))

print("LJ = %s" % cljff.energy( cljff.components().lj() ))

=================================

---------------------------------------------------------
Christopher Woods
+44 (0) 7786 264562
http://chryswoods.com

[Julien Michel]

Dear Sajid,

You can see an example script in

sire/corelib/test/SireIO/trajectory.py

This loads a pdb file that contains many water molecules. The first water molecule is manually assigned charges and LJ parameters, and the same parameters are assigned to every other water molecule.

For flexible molecules you need to setup intramolecular force field terms, and there is a lot more work to do to also set up connectivity, bonds, angles, dihedrals, excluded pairs.

You can see how this is done in C++ in sire/corelib/src/libs/SireIO/amber.* which contains code to parse amber topology files.

In general I recommend to generate a top/crd file,




--------------------------------------------------------------
Dr. Julien Michel,
Royal Society University Research Fellow
Room 263
School of Chemistry
University of Edinburgh
West Mains Road
Edinburgh, EH9 3JJ
United Kingdom
phone:   +44 (0)131 650 4797
http://www.julienmichel.net/
-------------------------------------------------------------

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

--
--------------------------------------------------------------
Dr. Julien Michel,
Royal Society University Research Fellow
Room 263
School of Chemistry
University of Edinburgh
West Mains Road
Edinburgh, EH9 3JJ
United Kingdom
phone:   +44 (0)131 650 4797
http://www.julienmichel.net/
-------------------------------------------------------------

[sajid iqbal]

Dear Dr. Christopher and Dr. Julien,

Thank you a lot. It works like a charm.

Regards,

Sajid

--
You received this message because you are subscribed to a topic in the Google Groups "Sire Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sire-users/q_x-ggmjqIo/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sire-users+...@googlegroups.com.

[sajid iqbal]

Dear Julien,

I agree with you regarding uploading amber topology and coordinate files. However, is there any possibility if we can modify the atoms in given molecule and then perform intra-molecule moves?

I would like to change the methyl group (present in the beginning of molecule) with hydrogen. Can I perform intra-moves without any need of defining parameters? I have attached the files. Please inform me if I can do it.

thank you

On Thu, Apr 21, 2016 at 4:19 AM, 'Julien Michel' via Sire Users <sire-...@googlegroups.com> wrote:

--
You received this message because you are subscribed to a topic in the Google Groups "Sire Users" group.
To unsubscribe from this topic, visit https://groups.google.com/d/topic/sire-users/q_x-ggmjqIo/unsubscribe.
To unsubscribe from this group and all its topics, send an email to sire-users+...@googlegroups.com.

[Julien Michel]

Dear Sajid,

It should be possible to run intra molecular move on flexible molecules without defining an internal force field and intramolecular non bonded force field. However the moves will generate distorted structures since changes in internal energy will not be evaluated.

There is to my knowledge currently no python script that sets up a full forcefield for flexible molecules. However if should be possible to do this by using the code in SireIO/amber.cpp as example to write a new python script. This will require some time and expertise in python but may be worth it if you do not wish to read a top/crd file.

Best wishes,

Julien

Sent from my iPhone

<box101.inpcrd>

<box101.prmtop>

[sajid iqbal]

Dear Julien,

Thank you a lot! I am able to change the the desired atoms and consequently generated the top/crd files before loading files on Sire. 

Regards,

Sajid