Having an optimized geometry of the ligands before MD simulations is required to have a good relative spatial arrangement of atoms in molecule, such as the right geometry of rings, 3D geometry of chiral centers, etc. However, a costly high performance geometry optimization may not be required, since this geometry optimization will be in the gas phase, which is different from the geometry of ligand in solution or in complex with a receptor.

Regarding the parametrization of the ligand before MD simulations, this depends on the force field you chose for simulations. You may chose doing parametrization using high level Quantum calculations using Gaussian, for example, however the potential form should be the same as the force field you chose for MD simulations, for example in your case amber. Alternatively, you may also chose Antechamber module with General Amber Force Field, which is designed to cover most pharmaceutical molecules and is compatible with the amber force fields in such a way that the two can be mixed during a simulation: it uses a simple harmonic function form for bonds and angles.The charge methods used can be HF/6-31G* RESP or AM1-BCC.

You can have a look at these three references:

J. Comp. Chem., 2004, 25, 1157.

J. Phys. Chem, 1993, 10269.

J. Comp. Chem., 2000, 21, 132.

Having an optimized geometry of the ligands before MD simulations is required to have a good relative spatial arrangement of atoms in molecule, such as the right geometry of rings, 3D geometry of chiral centers, etc. However, a costly high performance geometry optimization may not be required, since this geometry optimization will be in the gas phase, which is different from the geometry of ligand in solution or in complex with a receptor.

Regarding the parametrization of the ligand before MD simulations, this depends on the force field you chose for simulations. You may chose doing parametrization using high level Quantum calculations using Gaussian, for example, however the potential form should be the same as the force field you chose for MD simulations, for example in your case amber. Alternatively, you may also chose Antechamber module with General Amber Force Field, which is designed to cover most pharmaceutical molecules and is compatible with the amber force fields in such a way that the two can be mixed during a simulation: it uses a simple harmonic function form for bonds and angles.The charge methods used can be HF/6-31G* RESP or AM1-BCC.

You can have a look at these three references:

J. Comp. Chem., 2004, 25, 1157.

J. Phys. Chem, 1993, 10269.

J. Comp. Chem., 2000, 21, 132.

To add to Hiqmet: I agree that a high-level ligand optimization may not be needed, and, as long as they're rather simple organic molecules, semi-empirical methods shuold do just fine.

For ligand optimization before docking, it's even more complicated. There is a study (http://dx.doi.org/10.1021/jm030563w) that show that around 40% of ligands in the crystal structures have strain (the difference between the docked conformation and the solvent conformation) of more than 5 kcal/mol and 10% have strain of more than 9 kcal/mol. That shows you that ligands do not always bind to proteins in the same conformation, they're founds in the solution. If your docking algorithm can't handle ligand flexibility, there is a good chance you'll miss the correct pose and long minimizations of the ligand using some high-level QM does not give you any advantage.

You should always geometry-optimize your system with the same force field that you will then use for MD. If you run QM before FF-geometry optimization to get a low energy conformer and partial carges, that's fine.

With regard to the determination of partial atomic charges, optimization can be performed. But it is better to define them as an average over many conformations. In such cases I do a few short Hybrid MC simulations.

As to the MD, the optimization is needed only in the case if the original model is very rough. For example, if there are overlapping atoms. In this case it is sufficient to make several tens optimization steps. If the gradient is smaller that 1 the optimization is definitely not needed.

Moreover MD at low temperature is itself a good optimization method.

FF optimization gives a you kind of quality control measure, too. If your FF optimization drives you too far away from you QC conformer(s), then you might have assigned wrong atom types to your molecule.

Personally, I think the best way to handle ligand is to derive RESP charges (which in the same time optimizes its geometry) by RED Server, http://q4md-forcefieldtools.org/REDS/

Hello Nitin,

I've posted an answer to this question on the Amber listserve. However, for the sake of including here into this RG group discussion, I have copied-pasted it below. I think you will see common themes developing from everyone's answers.

---- cut-paste---

Hello Nitin,

Yes, in general one should optimize a ligand's geometry prior to an MD simulation. If the force-field parameters are present and trustworthy, then performing an MM minimization is sufficient. This is important so that you do not introduce artificially large forces, arising from the molecule residing at a high energy state on the MM's potential energy surface, at the start of your simulation.

However, if the parameters are not present (e.g. torsion terms, partial atomic charges) then you will need to optimize the ligand using QM or Antechamber (i.e. sqm). Most Amber force fields were developed with geometries optimized at HF/6-31G(d), but many use the Antechamber/sqm method for optimizing ligands, which I believe is semiempirical. Both methods are comparable, in the sense that they optimize the geometry, but differ in the theory used to do so. Different theories will lead to different optimized geometries, sometimes they have relatively minor differences (e.g. bond distances) and sometimes they are more significant (e.g. different rotamers).

Once you have the optimized conformation, then you use it as an input to determine partial atomic charges. There is quite a bit of literature for this already present for Amber, including using multiple conformations to determine a single charge set for a molecule.

It is reasonable/necessary to optimize a ligand's conformation that you obtained from a docking pose. If you used flexible docking methods, then note that this conformation may change due to the absence of the binding pocket.

As for websites/methods that you can use, Antechamber's can perform an optimization (via sqm) and to determine partial atomic charges. The R.E.D. website will also allow you to do both in a manner that is balanced with the Amber force fields. You can also have the option of downloading a number of QM software packages, such as Gamess and Gaussian. I would suggest doing some more reading regarding partial atomic charge determination and ligand/residue model building (e.g. Parm94, Gaff, Glycam06, R.E.D., Wolf2Pack, etc. papers).

Best regards,

Karl

Thanks for the replies

@HIQMET @BARTOZ as both you have mentioned that semi empirical methods are fine then it will be fine if i work with sqm which according to @KARL is a semi empirical method

for parametrization i can use AM1-BCC as it was parametrized to reproduce HF/6-31G* RESP charges and according to AMBER manual is recommended in large-scale calculations because of its efficiency

@BARTOZ as you mentioned that 40% of ligands in the crystal structures have strain of more than 5 kcal/mol and 10% have strain of more than 9 kcal/mol then it is a good approach to minimize ligand before docking. Moreover, when we optimise the conformation of the ligand we remove the strain and in my thinking that will deviate the result from the crystal structure. So in both cases, i.e. docking and MD after optimisation, comparison of the result with crystal pose doesn't make much sense.

@YUN SHI it will take me sometime to learn how to use RED server so for time being i will not be using it

@karl

After I run following command

antechamber -i ligand.pdb -fi pdb -o ligand.mol2 -fo mol2 -c bcc -s 2 -nc 0 I get sqm.pdb i.e. pdb file of the ligand generated after SQM, in addition to ligand.mol2.

1) the structure of sqm.pdb is different from original.pdb

2) if ligand.mol2 is close to original structure instead of sqm.pdb then what is the use of sqm.pdb file. Am I right in assuming that geometry optimization has not been achieved in real sense.

3) will it be a good idea to use sqm.pdb instead of ligand.mol2 to combine with receptor in order to form complex? i.e. UNK = loadpdb sqm.pdb instead of UNK = loadmol2 ligand.mol2

Also it will be really helpful if you can pass me some references related to partial atomic charge determination and ligand/residue model building.

Well it is hard to say because you are ultimately going to use a force-field to compute the conformations and dynamics of your system. I believe that sqm using the AM1 Hamiltonian for its calculations, which is typically pretty good for "normal" systems. You don't say how you obtained the original structure (i.e. build in a program, extracted from a crystal, etc.), so I can't say which structure I would preferentially use. (I am pretty sure that the -o ligand.mol2 is going to have the same coordinates as the -i ligand.pdb - you can verify this by looking into the files to see if they are the same.)

Ultimately, the MD simulation, IF given enough time, will provide you with a population of different conformations and configuration (interactions between ligand and receptor). The difficult thing is running the simulations long enough, which is extremely difficult to do for large systems and is dependent upon the force field and other approximations made within the theories. One thing that many researchers do is to run multiple long MD runs using the same and/or different input structures to show that their simulations converge. Statistically speaking the observable that you are interested should be converged to a value (e.g. interaction energy) or to several values (e.g. rotamer populations). For example, you could run two simulations, one with each of the ligand's conformations and then observe what happens. If the ligand is small and inflexible, you may observe that they converge to a similar geometry; if the ligand is large and flexible, you may see the simulations populate different geometries.

Ultimately, you have to validate your model building and methodology. One way to do this is the show that they reproduce experimental observables of your system. If you don't have experimental observables for your system, and you are using the MD theory to make predictions, then you should show that your methodology works in similar systems that have such data. I would suggest searching for similar studies in journals with high impact factors and see what they did.

Here are some references to get you started. This is by no means comprehensive, and is intended to give you a place to start your own literature search and review.

Best regards,

Karl

1. Bayly, C. I.; Cieplak, P.; Cornell, W. & Kollman, P. A. A well-behaved electrostatic potential based method using charge restraints for deriving atomic charges: the RESP model The Journal of Physical Chemistry, 1993, 97, 10269-10280

2. Application of RESP Charges to Calculate Conformational Energies, Hydrogen Bond Energies, and Free Energies of Solvation Journal of the American Chemical Society, 1993, 115, 9620-9631

3. Dupradeau, F.-Y.; Cézard, C.; Lelong, R.; Stanislawiak, É.; Pêcher, J.; Delepine, J. C. & Cieplak, P. R.E.DD.B.: A database for RESP and ESP atomic charges, and force field libraries Nucleic Acids Research, 2008, 36, D360-D367

4. Dupradeau, F.-Y.; Pigache, A.; Zaffran, T.; Savineau, C.; Lelong, R.; Grivel, N.; Lelong, D.; Rosanski, W. & Cieplak, P. The R.E.D. tools: advances in RESP and ESP charge derivation and force field library building Phys Chem Chem Phys, 2010, 12, 7821-7839

5. Lijiang Yang; Chun-hu Tan; Meng-Juei Hsieh; Junmei Wang; Yong Duan; Cieplak, P.; Caldwell, J.; Kollman, P. & Luo, R. New-generation Amber united-atom force field Journal of Physical Chemistry B, 2006, 110, 13166-76

6. Vanommeslaeghe, K. & MacKerell, A. D. Automation of the CHARMM General Force Field (CGenFF) I: Bond Perception and Atom Typing Journal of Chemical Information and Modeling, 2012, 52, 3144-3154

7. Wang, J.; Cieplak, P. & Kollman, P. A. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J. Comput. Chem., 2000, 21, 1049-1074

8. Kirschner, K. N.; Yongye, A. B.; Tschampel, S. M.; Gonzalez-Outeirino, J.; Daniels, C. R.; Foley, B. L. & Woods, R. J. GLYCAM06: A Generalizable Biomolecular Force Field. Carbohydrates Journal of Computational Chemistry, 2008, 29, 622-655

9. Reith, D. & Kirschner, K. N. A modern workflow for force-field development -- Bridging quantum mechanics and atomistic computational models Computer Physics Communications, 2011, 182, 2184 - 2191

10. Shim, J. & MacKerell, Jr., A. D. Computational ligand-based rational design: role of conformational sampling and force fields in model development Med. Chem. Commun., 2011, 2, 356-370

11. Kuhn, B.; Gerber, P.; Schulz-Gasch, T. & Stahl, M. Validation and use of the MM-PBSA approach for drug discovery. J. Med. Chem., 2005, 48, 4040-4048

12. Alonso, H.; Bliznyuk, A. A. & Gready, J. E. Combining docking and molecular dynamic simulations in drug design. Med. Res. Rev., 2006, 26, 531-568

13. Boehr, D. D.; Nussinov, R. & Wright, P. E. The role of dynamic conformational ensembles in biomolecular recognition. Nat Chem Biol, 2009, 5, 789-796

Thanks @karl

i extracted the structure from the crystal structure.

I have attached the images and coordinates file to make a clear picture of the problem

so can i run one simulation with ligand.mol2 and another with sqm.pdb

thanks for replying with so much patience

i am really grateful to you

There are a few things that you should consider for you structure, which are hard to tell from the images.

1) What is the overall charge of the molecule. You have both an amine that can be protonated and a carboxylic acid group that can be deprotonated. The charge state of you molecule will have a large impact on the binding. To help you determine this you can a) See what the crystal structure paper says, b) look at the binding site to see if an obvious ligand-amino acid pairing would benefit from (de-)protonation of the ligand, c) consult an organic chemist.

2) Aside from the double bond that the visualization program displays, are the other carbon atoms within the ring saturated or unsaturated carbons? Make sure you have this right, and don't rely on what the visualization program/antechamber decides it is.

3) The ether side group is the one that should the largest conformational difference between the crystal and the sqm structure. Structures a and c are the same (as mentioned before), but differ only in the added hydrogens. My guess is that the MD simulations will sort out which side-chain conformation (both dihedral atoms defined by OS-CT-CT-CT) will be preferred, and will depend partly on the binding site geometry. Running two simulations, using the two structure as input, would be the safe thing to do if you can afford it.

Remember that the trustworthiness of the conformation coming from the crystal structure is dependent upon the resolution of the crystal structure. It also represents a solid state solution (and typically at lower temperature than the body), while you will be doing a solution-state simulation. It is very reasonable to see some local conformational changes to occur in your solution-state, 310 K dynamics.

Best regards,

Karl

If you design the ligand from scratch all levels of geometry optimization (QM,MM) are useful. If you have already a crystal structure, geometry optimization can be done with the force field you use in MD simulation.

Anyway, please pay attention to what Karl Kirschner noted.

Thanks a lot @karl

@georgious thanks for the reply

does anyone know how to use the Gromac to do the geometry optimization of peptide and the energy decomposition analysis?

Thank you very much!