[gmx-users] Cyclohexane simulation problem
Justin Lemkul
jalemkul at vt.edu
Sun May 31 17:08:24 CEST 2015
On 5/31/15 6:00 AM, mohsen shahlaei wrote:
> Thank you Justin for prompt reply.I have been used coulombtype = Cut-off in my mdp file for Cyclohexane energy minimization but the results sound odd.I used cyclohexane density (0.824 g/ml) for calculating number of CHX molecules in a 10 nm3 (about 56 molecules). Then I used following command to make a box.gmx insert-molecules -ci chx.gro -nmol 56 -box 2.154434 2.154434 2.154434 -o chx_box.grochx.gro and chx.itp that I used are:
> =============chx.itp ==================;
> ;
> ; This file was generated by PRODRG version AA100323.0717
> ; PRODRG written/copyrighted by Daan van Aalten
> ; and Alexander Schuettelkopf
> ;
> ; Questions/comments to dava at davapc1.bioch.dundee.ac.uk
> ;
> ; When using this software in a publication, cite:
> ; A. W. Schuettelkopf and D. M. F. van Aalten (2004).
> ; PRODRG - a tool for high-throughput crystallography
> ; of protein-ligand complexes.
> ; Acta Crystallogr. D60, 1355--1363.
> ;
> ;
>
> [ moleculetype ]
> ; Name nrexcl
> chx 3
>
> [ atoms ]
> ; nr type resnr resid atom cgnr charge mass
> 1 CH1 1 chx CAA 1 0.000 14.0270
> 2 CH1 1 chx CAB 1 0.000 14.0270
> 3 CH1 1 chx CAE 1 0.000 14.0270
> 4 CH1 1 chx CAF 1 0.000 14.0270
> 5 CH1 1 chx CAD 1 0.000 14.0270
> 6 CH1 1 chx CAC 1 0.000 14.0270
>
> [ bonds ]
> ; ai aj fu c0, c1, ...
> 1 2 2 0.152 5430000.0 0.152 5430000.0 ; CAA CAB
> 1 6 2 0.152 5430000.0 0.152 5430000.0 ; CAA CAC
> 2 3 2 0.152 5430000.0 0.152 5430000.0 ; CAB CAE
> 3 4 2 0.152 5430000.0 0.152 5430000.0 ; CAE CAF
> 4 5 2 0.152 5430000.0 0.152 5430000.0 ; CAF CAD
> 5 6 2 0.152 5430000.0 0.152 5430000.0 ; CAD CAC
>
> [ pairs ]
> ; ai aj fu c0, c1, ...
> 1 4 1 ; CAA CAF
> 2 5 1 ; CAB CAD
> 3 6 1 ; CAE CAC
>
> [ angles ]
> ; ai aj ak fu c0, c1, ...
> 2 1 6 2 109.5 520.0 109.5 520.0 ; CAB CAA CAC
> 1 2 3 2 109.5 520.0 109.5 520.0 ; CAA CAB CAE
> 2 3 4 2 109.5 520.0 109.5 520.0 ; CAB CAE CAF
> 3 4 5 2 109.5 520.0 109.5 520.0 ; CAE CAF CAD
> 4 5 6 2 109.5 520.0 109.5 520.0 ; CAF CAD CAC
> 1 6 5 2 109.5 520.0 109.5 520.0 ; CAA CAC CAD
>
> [ dihedrals ]
> ; ai aj ak al fu c0, c1, m, ...
> 3 2 1 6 1 0.0 5.9 3 0.0 5.9 3 ; dih CAE CAB CAA CAC
> 5 6 1 2 1 0.0 5.9 3 0.0 5.9 3 ; dih CAD CAC CAA CAB
> 4 3 2 1 1 0.0 5.9 3 0.0 5.9 3 ; dih CAF CAE CAB CAA
> 5 4 3 2 1 0.0 5.9 3 0.0 5.9 3 ; dih CAD CAF CAE CAB
> 6 5 4 3 1 0.0 5.9 3 0.0 5.9 3 ; dih CAC CAD CAF CAE
> 1 6 5 4 1 0.0 5.9 3 0.0 5.9 3 ; dih CAA CAC CAD CAF
> =================================chx.gro==================== PRODRG COORDS
> 6
> 1chx C AA 1 0.155 0.105 0.177
> 1chx C AB 2 0.100 0.247 0.184
> 1chx C AE 3 0.182 0.344 0.100
> 1chx C AF 4 0.330 0.339 0.137
> 1chx C AD 5 0.385 0.196 0.130
> 1chx C AC 6 0.303 0.100 0.214
> 0.48500 0.44400 0.31400===========================================================
As I suspected before, you have a malformed coordinate file. Fix it before
doing anything.
http://manual.gromacs.org/online/gro.html
> Then I minimized the energy of generated box as follows:gmx grompp -f minim.mdp -c chx_box.gro -p chx_box.top -o em.tpr -maxwarn 1gmx mdrun -deffnm em -v
>
>
> The warning I got in grompp runing was
> WARNING 1 [file chx_box.top, line 58]:
> 336 non-matching atom names
> atom names from chx_box.top will be used
> atom names from chx_box.gro will be ignored
>
>
> Removing all charge groups because cutoff-scheme=Verlet
> Analysing residue names:
> There are: 56 Other residues
> Analysing residues not classified as Protein/DNA/RNA/Water and splitting into groups...
> Number of degrees of freedom in T-Coupling group rest is 1005.00
> This run will generate roughly 2 Mb of data
>
> There was 1 note
>
> There was 1 warning
>
> -------------------------------------------------------
> Program gmx, VERSION 5.0.4
> Source code file: /home/mohsen/Downloads/gromacs-5.0.4/src/gromacs/gmxpreprocess/grompp.c, line: 2087
>
> Fatal error:
> Too many warnings (1), gmx terminated.
>
> So I used -maxwarn 1
Generally dangerous. Perhaps not relevant here, but you shouldn't get into this
habit. Fix the problems in the input rather than overriding warnings.
> The chx_box ITP is==========================chx_box.itp=====================;
> ; File 'topol.top' was generated
> ; By user: mohsen (1000)
> ; On host: mohsen-All-Series
> ; At date: Wed May 20 12:09:46 2015
> ;
> ; This is a standalone topology file
> ;
> ; Created by:
> ; GROMACS: gmx pdb2gmx, VERSION 5.0.4
> ; Executable: /usr/local/gromacs/bin/gmx
> ; Library dir: /usr/local/gromacs/share/gromacs/top
> ; Command line:
> ; gmx pdb2gmx -f pep.gro -o peptide.gro -water spce
> ; Force field was read from the standard Gromacs share directory.
> ;
> nrexcl=3
> ; Include forcefield parameters
> #include "gromos53a6.ff/forcefield.itp"
>
>
> ; Include Position restraint file
> #ifdef POSRES
> #include "posre.itp"
> #endif
>
> ; octanol position restraints
> #ifdef POSRES_CHX
> #include "system.itp"
> #endif
>
> ; Include water topology
> #include "gromos53a6.ff/spce.itp"
> #include "chx.itp"
>
>
> #ifdef POSRES_WATER
> ; Position restraint for each water oxygen
> [ position_restraints ]
> ; i funct fcx fcy fcz
> 1 1 1000 1000 1000
> #endif
>
> ; Include topology for ions
> #include "gromos53a6.ff/ions.itp"
>
> [ system ]
> ; Name
> Gromacs Runs On Most of All Computer Systems in
>
> [ molecules ]
> ; Compound #mols
> chx 56
> and the minim.mdp is as follows:===============================minim.itp=======================; minim.mdp - used as input into grompp to generate em.tpr
> integrator = steep ; Algorithm (steep = steepest descent minimization)
> emtol = 100.0 ; Stop minimization when the maximum force < 1000.0 kJ/mol/nm
> emstep = 0.01 ; Energy step size
> nsteps = 500000 ; Maximum number of (minimization) steps to perform
>
> ; Parameters describing how to find the neighbors of each atom and how to calculate the interactions
> 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)
> coulombtype = Cut-off ; Treatment of long range electrostatic interactions
> rcoulomb = 1 ; Short-range electrostatic cut-off
> rvdw = 1 ; Short-range Van der Waals cut-off
> pbc = xyz ; Periodic Boundary Conditions (yes/no)===========================================================
> the results of energy minimization is as follows:writing lowest energy coordinates.
>
> Back Off! I just backed up em.gro to ./#em.gro.19#
>
> Steepest Descents converged to Fmax < 100 in 1025 steps
> Potential Energy = 1.9984604e+03
> Maximum force = 9.3209557e+01 on atom 181
> Norm of force = 2.5035557e+01
>
> NOTE: 19 % of the run time was spent in pair search,
> you might want to increase nstlist (this has no effect on accuracy)
> As it can be seen the potential energy is positive. Please let me know why and solution please.
>
So? The energy of the system depends entirely upon the summation of all the
interactions. You have weakly interacting species - vdW only, no charges.
Bonded interactions are strictly positive. So if the sum of Ebonded + EvdW > 0
this is what you'll get. There is no reason to think this is inherently wrong.
A little uncommon for a condensed-phase system, and not something you would
see in an aqueous system, but for something like this it's probably fine.
-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
jalemkul at outerbanks.umaryland.edu | (410) 706-7441
http://mackerell.umaryland.edu/~jalemkul
==================================================
More information about the gromacs.org_gmx-users
mailing list