[gmx-users] Diagnosing + system blowing up
Mark Abraham
Mark.Abraham at anu.edu.au
Tue Jul 31 02:10:15 CEST 2012
On 31/07/2012 7:16 AM, Shima Arasteh wrote:
> Wow! That was a fabulous explanation given to me. Thanks dear Mark! :-)
>
> I regenerated the topol.top to check the correctness of input and FF files. Some output has been changed however a little bit. Now, I bring you all was needed again:
>
>> 1. rtp entry
>>
>> [ FVAL ]
> [ atoms ]
> CN C 0.357 0
> ON O -0.51 1
> H1 HA 0.100 2
> N NH1 -0.423 3
> HN H 0.333 4
> CA CT1 0.034 5
> HA HB 0.09 6
> CB CT1 -0.093 7
> HB HA 0.09 8
> CG1 CT3 -0.268 9
> HG11 HA 0.09 10
> HG12 HA 0.09 11
> HG13 HA 0.09 12
> CG2 CT3 -0.268 13
> HG21 HA 0.09 14
> HG22 HA 0.09 15
> HG23 HA 0.09 16
> C C 0.528 17
> O O -0.510 18
> [ bonds ]
> CN H1
> CN ON
> CN N
> N HN
> CA N
> CA HA
> CA C
> C O
> CA CB
> CB HB
> CB CG1
> CB CG2
> CG2 HG21
> CG2 HG22
> CG2 HG23
> CG1 HG11
> CG1 HG12
> CG1 HG13
>
> [ impropers ]
> CN N ON H1
>
>> 2. hdb entry
> In which file?
> *aminoacids.hdb
>
>
>> FVAL 6
>> 1 1 H1 CN N ON
> This should be generating H1 bonded to CN...
>
>> 1 1 HN N C CA
>> 1 5 HA CA N C CB
>> 1 5 HB CB CA CG1 CG2
>> 3 4 HG1 CG1 CB CA
>> 3 4 HG2 CG2 CB CA
>>
>> 3. N-terminal fragment
>> HETATM 1 CN FVAL 1 -0.721 1.600 1.249
>> HETATM 2 ON FVAL 1 -0.839 2.806 1.453
>> ATOM 3 N FVAL 1 -1.227 0.728 2.125
>> ATOM 4 CA FVAL 1 -1.918 1.159 3.323
>> ATOM 5 C FVAL 1 -1.969 2.678 3.410
>> ATOM 6 O FVAL 1 -0.931 3.335 3.447
>> ATOM 7 CB FVAL 1 -1.219 0.644 4.576
>> ATOM 8 CG1 FVAL 1 0.208 1.178 4.618
>> ATOM 9 CG2 FVAL 1 -1.976 1.118 5.812
>>
>> 4. pdb2gmx command
>> #pdb2gmx -f monomer.pdb -o monomer.gro -water tip3p -ter
>>
>> 5. pdb2gmx output
> Using the Charmm36-modified force field in directory ./charmm36-modified.ff
>
> Opening force field file ./charmm36-modified.ff/aminoacids.r2b
> Opening force field file ./charmm36-modified.ff/rna.r2b
> Reading monomer.pdb...
> Read 177 atoms
> Analyzing pdb file
> Splitting chemical chains based on TER records or chain id changing.
> There are 1 chains and 0 blocks of water and 24 residues with 177 atoms
>
> chain #res #atoms
> 1 ' ' 24 177
>
> All occupancy fields zero. This is probably not an X-Ray structure
> Opening force field file ./charmm36-modified.ff/atomtypes.atp
> Atomtype 1
> Reading residue database... (charmm36-modified)
> Opening force field file ./charmm36-modified.ff/aminoacids.rtp
> Residue 42
> Sorting it all out...
> Opening force field file ./charmm36-modified.ff/dna.rtp
> Residue 46
> Sorting it all out...
> Opening force field file ./charmm36-modified.ff/lipids.rtp
> Residue 82
> Sorting it all out...
> Opening force field file ./charmm36-modified.ff/rna.rtp
> Residue 86
> Sorting it all out...
> Opening force field file ./charmm36-modified.ff/aminoacids.hdb
> Opening force field file ./charmm36-modified.ff/dna.hdb
> Opening force field file ./charmm36-modified.ff/lipids.hdb
> Opening force field file ./charmm36-modified.ff/rna.hdb
> Opening force field file ./charmm36-modified.ff/aminoacids.n.tdb
> Opening force field file ./charmm36-modified.ff/dna.n.tdb
> Opening force field file ./charmm36-modified.ff/rna.n.tdb
> Opening force field file ./charmm36-modified.ff/aminoacids.c.tdb
> Opening force field file ./charmm36-modified.ff/dna.c.tdb
> Opening force field file ./charmm36-modified.ff/rna.c.tdb
>
> Back Off! I just backed up topol.top to ./#topol.top.2#
> Processing chain 1 (177 atoms, 24 residues)
> Identified residue FVAL1 as a starting terminus.
> Identified residue GLY24 as a ending terminus.
> 8 out of 8 lines of specbond.dat converted successfully
> Select start terminus type for FVAL-1
> 0: NH3+
> 1: NH2
> 2: None
> 2
> Start terminus FVAL-1: None
> Select end terminus type for GLY-24
> 0: COO-
> 1: COOH
> 2: CT2
> 3: CT3
> 4: None
> 0
> End terminus GLY-24: COO-
> Opening force field file ./charmm36-modified.ff/aminoacids.arn
> Opening force field file ./charmm36-modified.ff/dna.arn
> Opening force field file ./charmm36-modified.ff/rna.arn
> Checking for duplicate atoms....
> Now there are 24 residues with 360 atoms
> Making bonds...
> Number of bonds was 362, now 362
> Generating angles, dihedrals and pairs...
> Before cleaning: 925 pairs
> Before cleaning: 930 dihedrals
> Keeping all generated dihedrals
> Making cmap torsions...There are 22 cmap torsion pairs
> There are 930 dihedrals, 49 impropers, 644 angles
> 916 pairs, 362 bonds and 0 virtual sites
> Total mass 2510.906 a.m.u.
> Total charge 1.000 e
> Writing topology
>
> Back Off! I just backed up posre.itp to ./#posre.itp.1#
>
> Writing coordinate file...
>
> Back Off! I just backed up monomer.gro to ./#monomer.gro.1#
> --------- PLEASE NOTE ------------
> You have successfully generated a topology from: monomer.pdb.
> The Charmm36-modified force field and the tip3p water model are used.
> --------- ETON ESAELP ------------
>
>
>> 6. N-terminal fragment of the grompp -c input
> 1FVAL CN 1 3.039 1.904 2.416
> 1FVAL H1 2 3.086 1.871 2.334
This H1 atom is now in the right position for generation from the
correct .hdb, and pdb2gmx is not complaining of the long bond to H1,
which I think means you are using the correct .hdb now, and were not
doing so earlier. So my suggestion of a bug was unfounded.
> 1FVAL ON 3 3.027 2.025 2.436
> 1FVAL N 4 2.988 1.817 2.504
> 1FVAL HN 5 3.029 1.759 2.433
> 1FVAL CA 6 2.919 1.860 2.623
> 1FVAL HA 7 2.828 1.821 2.614
> 1FVAL CB 8 2.989 1.808 2.749
> 1FVAL HB 9 2.991 1.708 2.747
> 1FVAL CG1 10 3.132 1.862 2.753
> 1FVAL HG11 11 3.178 1.828 2.835
> 1FVAL HG12 12 3.181 1.831 2.672
> 1FVAL HG13 13 3.130 1.962 2.754
> 1FVAL CG2 14 2.913 1.856 2.872
> 1FVAL HG21 15 2.959 1.822 2.954
> 1FVAL HG22 16 2.911 1.956 2.874
> 1FVAL HG23 17 2.820 1.821 2.869
> 1FVAL C 18 2.914 2.012 2.632
> 1FVAL O 19 3.018 2.077 2.636
>
> ... and this shows atom H1 around 0.1nm from atom C. I can think of no
> legitimate reason how this could occur, given your stated .hdb, but I
> seem to recall an old .hdb version you showed had a line
>
> 1 1 H1 C N ON
>
> rather than the correct
>
> 1 1 H1 CN N ON
>
> ********* I checked it, there was the correct one.
>
> and this would explain the weird H1 position perfectly (but not how it still gets bonded to atom 1). I'd like you to double check that ./charmm36-modified.ff/aminoacids.hdb has (only) the correct FVAL entry. Also, I'd like to see the first 30 or so lines of the [bonds] section of this [moleculetype] in topol.top, so we can really see what bonds are generated.
>
> * First 50 lines of the [bonds] section:
These are now normal. I suspect they were not earlier.
>
> [ bonds ]
> ; ai aj funct c0 c1 c2 c3
> 1 2 1
> 1 3 1
> 1 4 1
> 4 5 1
> 4 6 1
> 6 7 1
> 6 8 1
> 6 18 1
> 8 9 1
> 8 10 1
> 8 14 1
> 10 11 1
> 10 12 1
> 10 13 1
> 14 15 1
> 14 16 1
> 14 17 1
> 18 19 1
> 20 21 1
> 20 22 1
> 22 23 1
> 22 24 1
> 22 29 1
> 24 25 1
> 24 26 1
> 24 27 1
> 27 28 1
> 29 30 1
> 29 31 1
> 31 32 1
> 31 33 1
> 33 34 1
> 33 35 1
> 33 48 1
> 35 36 1
> 35 37 1
> 35 38 1
> 38 39 1
> 38 40 1
> 38 44 1
> 40 41 1
> 40 42 1
> 40 43 1
> 44 45 1
> 44 46 1
> 44 47 1
> 48 49 1
> 48 50 1
> 50 51 1
> 50 52 1
>
>
> I think you may have encountered a new bug in pdb2gmx.
>
>> 1FVAL O 19 3.019 2.076 2.635
>>
>> 7.grompp command line
>> # grompp -f ions.mdp -c monomer_solv.gro -p topol.top -o ions.tpr
>>
>>
>> 8. grompp output
>> Generated 21528 of the 21528 non-bonded parameter combinations
>> Generating 1-4 interactions: fudge = 1
>> Generated 18355 of the 21528 1-4 parameter combinations
>>
>> ERROR 1 [file topol.top, line 414]:
>> No default Bond types
>>
>>
>> ERROR 2 [file topol.top, line 1698]:
>> No default U-B types
>>
>>
>> ERROR 3 [file topol.top, line 1699]:
>> No default U-B types
>>
>>
>> ERROR 4 [file topol.top, line 2345]:
>> No default Proper Dih. types
>>
>>
>> ERROR 5 [file topol.top, line 2346]:
>> No default Proper Dih. types
>>
>> ERROR 6 [file topol.top, line 3278]:
>> No default Improper Dih. types
> All this probably follows on from the confused state of bonding to H1,
> e.g. if CN-H1 is defined from the .rtp bond, and there's now a C-H1 from
> the possibly erroneous .hdb construction. What is line 414 of topol.top?
>
> * line 414 of topol.tp:
> [bonds]
> 1 2 1
That's what it is in the new version, not the old version that actually
had the errors.
What does grompp have to say now?
Mark
More information about the gromacs.org_gmx-users
mailing list