On 2/13/15 6:03 AM, Jennifer Vo wrote:
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.
You'll have to be much more specific than "goes wrong with the ligands." What's
happening? Does EM fail or is the output just somehow screwed up? The ligand
files are both set up within different boxes, so perhaps the construction of
your system is incorrect somehow, but again, without specifics, I'm just guessing.
-Justin
--
==================================================
Justin A. Lemkul, Ph.D.
Ruth L. Kirschstein NRSA Postdoctoral Fellow
Department of Pharmaceutical Sciences
School of Pharmacy
Health Sciences Facility II, Room 629
University of Maryland, Baltimore
20 Penn St.
Baltimore, MD 21201
jalem...@outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul
==================================================
--
Gromacs Users mailing list
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a
mail to gmx-users-requ...@gromacs.org.
