[gmx-users] multiple ligands topology

[Jennifer Vo]

Dear All,
I am running a simulation for a systems including two chains of proteins
and two ligands using amber 99SB ff.
My topol.top is

#include "amber99sb.ff/forcefield.itp"
#include "my_ligand_atomtypes.itp"

; Include chain topologies
#include "A.itp"
#ifdef POSRES
#include "posre_A.itp"
#endif
#include "B.itp"
#ifdef POSRES
#include "posre_B.itp"
#endif

; Include custom ligand topologies
#include "npd.itp"
#ifdef POSRES_LIG
#include "posre_npd.itp"
#endif
; Include custom ligand topologies
#include "ir3_em.itp"
; Include water topology
#include "amber99sb.ff/tip3p.itp"

#ifdef POSRES_WATER
; Position restraint for each water oxygen
[ position_restraints ]
;  i funct       fcx        fcy        fcz
   1    1       1000       1000       1000
#endif
; Include generic ion topology
#include "amber99sb.ff/ions.itp"

[ system ]
; Name
BcSIRED in water
[ molecules ]
; Compound        #mols
A             1
B             1
NPD                   1
IR3                   9
SOL         123406
NA               24

then I have my_ligand_atomtypes.itp
[ atomtypes ]
;name   bond_type     mass     charge   ptype   sigma         epsilon
Amb
 CA       CA          0.00000  0.00000   A     3.39967e-01   3.59824e-01 ;
1.91  0.0860
 H4       H4          0.00000  0.00000   A     2.51055e-01   6.27600e-02 ;
1.41  0.0150
 HA       HA          0.00000  0.00000   A     2.59964e-01   6.27600e-02 ;
1.46  0.0150
 C        C           0.00000  0.00000   A     3.39967e-01   3.59824e-01 ;
1.91  0.0860
 O        O           0.00000  0.00000   A     2.95992e-01   8.78640e-01 ;
1.66  0.2100
 N        N           0.00000  0.00000   A     3.25000e-01   7.11280e-01 ;
1.82  0.1700
 H        H           0.00000  0.00000   A     1.06908e-01   6.56888e-02 ;
0.60  0.0157
 N*       N*          0.00000  0.00000   A     3.25000e-01   7.11280e-01 ;
1.82  0.1700
 CT       CT          0.00000  0.00000   A     3.39967e-01   4.57730e-01 ;
1.91  0.1094
 H2       H2          0.00000  0.00000   A     2.29317e-01   6.56888e-02 ;
1.29  0.0157
 H1       H1          0.00000  0.00000   A     2.47135e-01   6.56888e-02 ;
1.39  0.0157
 OH       OH          0.00000  0.00000   A     3.06647e-01   8.80314e-01 ;
1.72  0.2104
 HO       HO          0.00000  0.00000   A     0.00000e+00   0.00000e+00 ;
0.00  0.0000
 OS       OS          0.00000  0.00000   A     3.00001e-01   7.11280e-01 ;
1.68  0.1700
 P        P           0.00000  0.00000   A     3.74177e-01   8.36800e-01 ;
2.10  0.2000
 O2       O2          0.00000  0.00000   A     2.95992e-01   8.78640e-01 ;
1.66  0.2100
 CK       CK          0.00000  0.00000   A     3.39967e-01   3.59824e-01 ;
1.91  0.0860
 H5       H5          0.00000  0.00000   A     2.42146e-01   6.27600e-02 ;
1.36  0.0150
 NB       NB          0.00000  0.00000   A     3.25000e-01   7.11280e-01 ;
1.82  0.1700
 CB       CB          0.00000  0.00000   A     3.39967e-01   3.59824e-01 ;
1.91  0.0860
 N2       N2          0.00000  0.00000   A     3.25000e-01   7.11280e-01 ;
1.82  0.1700
 NC       NC          0.00000  0.00000   A     3.25000e-01   7.11280e-01 ;
1.82  0.1700
 CQ       CQ          0.00000  0.00000   A     3.39967e-01   3.59824e-01 ;
1.91  0.0860
 ; IR3_GMX.top created by acpype (Rev: 403) on Wed Feb  4 11:55:46 2015
; Include forcefield parameters
[ atomtypes ]
 ;name   bond_type     mass     charge   ptype   sigma
epsilon       Amb
 C        C           0.00000  0.00000   A     3.39967e-01   3.59824e-01 ;
1.91  0.0860
 H        H           0.00000  0.00000   A     1.06908e-01   6.56888e-02 ;
0.60  0.0157
 N        N           0.00000  0.00000   A     3.25000e-01   7.11280e-01 ;
1.82  0.1700

My first ligand gro file
GRoups of Organic Molecules in ACtion for Science
   73
    1NPD   OA22    1  -0.282   0.251  -5.719
    1NPD   OA23    2  -0.142   0.366  -5.542
    1NPD   OA24    3  -0.039   0.314  -5.771
    1NPD   P'A2    4  -0.145   0.285  -5.669
    1NPD     PA    5  -0.033  -0.056  -5.010
    1NPD     PN    6  -0.130  -0.245  -4.804
    1NPD    O3P    7  -0.126  -0.116  -4.893
    1NPD    N1A    8  -0.692   0.110  -5.840
    1NPD    N1N    9  -0.539  -0.355  -4.443
    1NPD    C2A   10  -0.613   0.000  -5.836
...
    1NPD    H8A   73  -0.376   0.292  -5.440
   1.94588     0.88647     1.08377

My second ligand gro file

IR3_GMX.gro created by acpype (Rev: 403) on Wed Feb  4 11:55:46 2015
22
    1IR3    CAA    1  10.182  10.576  10.387
    1IR3    HAB    2  10.211  10.640  10.305
    1IR3    HAA    3  10.075  10.577  10.396
....
   1IR3    HAF   22  10.581  10.138  10.385
  20.72517  20.72517  20.72517

but the system can`t go under energy minimization. This is my em.mdp

title           = Minimization  ; Title of run
integrator      = steep         ; Algorithm (steep = steepest descent
minimization)
emtol           = 1000.0        ; Stop minimization when the maximum force
< 1000.0 kJ/mol
emstep          = 0.01      ; Energy step size
nsteps          = 500000                ; Maximum number of (minimization)
steps to perform
energygrps      = Protein NPD IR3       ; Which energy group(s) to write to
disk
nstlist             = 1             ; Frequency to update the neighbor list
and long range forces
cutoff-scheme       = Verlet
ns_type             = grid              ; Method to determine neighbor list
(simple, grid)
rlist               = 1.0               ; Cut-off for making neighbor list
(short range forces)
coulombtype         = PME               ; Treatment of long range
electrostatic interactions
rcoulomb            = 1.0               ; long range electrostatic cut-off
rvdw                = 1.0               ; long range Van der Waals cut-off
pbc                     = xyz           ; Periodic Boundary Conditions

I change emtol from 1000, then 500, then 100. No error generated but the
pdb output file goes wrong with the ligands. If I do the same procedure
with seperated ligand (only one [ atomtypes] in my_ligand_atomtypes.itp),
things goes well . Could you please help me where is the problem?
Many thanks in advance.
Jennifer