[gmx-users] md with multiple ligands

[Justin Lemkul]

On 2/18/13 10:21 AM, Ansuman Biswas wrote:
> ---------------------------- Original Message ----------------------------
> Subject: md with multiple ligands
> From:    ansuman at localhost
> Date:    Mon, February 18, 2013 8:28 pm
> To:      gmx-users at gromacs.org
> --------------------------------------------------------------------------
>
> Message: 8
>>> Date: Wed, 13 Feb 2013 08:28:00 -0500
>>> From: Justin Lemkul <[hidden email]>
>>> Subject: Re: [gmx-users] md with multiple ligands
>>> To: Discussion list for GROMACS users <[hidden email]>
>>> Message-ID: <[hidden email]>
>>> Content-Type: text/plain; charset=ISO-8859-1; format=flowed
>>>
>>>
>>>
>>> On 2/13/13 5:24 AM, Ansuman Biswas wrote:
>>>>
>>>> Dear gromacs users,
>>>>
>>>> I wish to run a protein ligand md in Gromacs 4.5.3 with charmm force
>>>> field. Its a dimeric protein with a MG ion and 2 ligands bound to two
>>>> active sites. I wish to run the dimer and there would be 4 ligands. I
>>>> generated the topology files (.itp) of ligands using SWISSPARAM. Before
>>>> running the energy minimization I changed the topol.top file accordingly
>>>> as mentioned in the SWISSPARAM website. I added all the itp files
>>>> together.But while running GROMPP I got an error message like,
>>>>
>>>>
>>>> WARNING 1 [file ADP.itp, line 24]:
>>>>      Overriding atomtype NPYL
>>>>
>>>>
>>>> WARNING 2 [file ADP.itp, line 25]:
>>>>      Overriding atomtype C5A
>>>>
>>>>
>>>> WARNING 3 [file ADP.itp, line 26]:
>>>>      Overriding atomtype N5B
>>>>
>>>>
>>>> WARNING 4 [file ADP.itp, line 27]:
>>>>      Overriding atomtype C5B
>>>>
>>>>
>>>> Generated 25425 of the 25425 non-bonded parameter combinations
>>>> Generating 1-4 interactions: fudge = 1
>>>> Generated 22285 of the 25425 1-4 parameter combinations
>>>>
>>>> -------------------------------------------------------
>>>> Program grompp, VERSION 4.5.3
>>>> Source code file: toppush.c, line: 1071
>>>>
>>>> Fatal error:
>>>> Atoms in the .top are not numbered consecutively from 1 (rather, atomnr
>>>> =
>>>> 1, while at->nr = 40)
>>>>
>>>> I know that the problem is related to the merging of 4 itp files.
>>>> Please let me how I should do that.
>>>>
>>>
>>> The ligands probably use common atom types and each .itp file likely has
>>> an
>>> [atomtypes] directive that introduces these new types.  Probably the
>>> easiest
>>> thing to do is create a merged [atomtypes] directive that lists all of the
>>> necessary atom types so you don't list them separately in each file.  One
>>> approach that I have used in the past is to create an .itp file that only
>>> has
>>> atom types in it, then call something like:
>>>
>>> #include "charmm27.ff/forcefield.itp"
>>>
>>> ; only has [atomtypes]
>>> #include "my_ligand_atomtypes.itp"
>>>
>>> ; these don't introduce [atomtypes] any more
>>> #include "ligand1.itp"
>>> #include "ligand2.itp"
>>>
>>> -Justin
>> But according to SWISSPARAM website I can only add one ligand.itp file.
> «  [hide part of quote]
>
> Just a guess, but it's probably because each ligand.itp file introduces
> new atom
> types which give the conflicts above.
>
>> So, is it possible to add ligand1 and ligand2 itp files?
>>
>
> Sure, my method should work.  Please try it.  I have done the exact same
> thing
> with Amber topologies from acpype.
>
>> Will it be ok if I merge all the ligand pdbs and then get the combined itp
>> file using that from SWISSPARAM and follow the protocol for one ligand?
>>
>
> Merging the coordinate files will probably not work.  If you have multiple
> molecules, you need multiple topologies.
>
> Justin
>
>
> Well, I have tried using the method you mentioned; i.e using 2 separate
> itp files and use one common atomtype at the beginning. But it seemed that
> gromacs was unable to read the next lig1 and lig2 itp files without
> atomtypes option. It gave error message. So I tried merging two itp files

Well, what was the error?  If you want free help, you have to make it easy to 
help you.

> from swissparam so that only unique parameters remain in the combined.itp
> file. I could run Energy Minimization only with steepest descent but not
> with CG as CG could not reduce the maximum force beyond 5th order. After

What is 5th order?  You mean you have an Fmax > 10^5?  That would indicate 
significant instability.

> that when I tried running NVT, I found Lincs error and saw several bonds
> rotated beyond 30 deg, as well as domain decomposition error:
>          constraint lengths are too long compared to the domain
> decomposition cell size
>

If my assumption about the Fmax > 10^5 is correct, then this is not unexpected.

> *******More on that when I opened the em.pdb (em output) in pymol I found
> all the ligand atoms were jumbled up together.
>

This outcome indicates significant problems with the ligand topology(ies). 
Without the content of the ligand topologies, it's impossible to help you 
resolve this.

-Justin

-- 
========================================

Justin A. Lemkul, Ph.D.
Research Scientist
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin

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