Re: [gmx-users] doubt in mdrun
Hi Shine, The amount of data is due to a very long simulation (1 microsecond!) coupled to a very frequent dumping of snapshots (every 1000 frames). You should question whether you need such frequent storage of the trajectory (e.g., for analysis of diffusion) or if less frequent will do just fine. Also, question whether you need such a long simulation or a shorter one would be sufficient. In any case, you could do a shorter run, analyze the trajectory and continue the simulation as described in http://www.gromacs.org/Documentation/How-tos/Extending_Simulations . Cheers, Pablo Dr. Pablo Englebienne Postdoctoral Researcher *TU Delft / 3mE / Process Energy* /Engineering Thermodynamics (ETh) group/ Building 46 Leeghwaterstraat 44, room 030 2628 CA Delft The Netherlands *T* +31 (0)15 27 86662 tel:+31152786662 *E* p.englebie...@tudelft.nl mailto:p.englebie...@tudelft.nl On 05/01/2013 07:29 AM, Shine A wrote: Sir, I studying the dynamics of a peptide in explicit solvent model.But during the mdrun I got the message like this. NOTE 1 [file md.mdp]: This run will generate roughly 1177130 Mb of data why the run generating this much amount of data? The md.mdp file I used is shown below title = OPLS Lysozyme MD ; Run parameters integrator = md; leap-frog integrator nsteps = 5 ; 2 * 5 = 100 ps, 1000 ns dt = 0.002 ; 2 fs ; Output control nstxout = 1000 ; save coordinates every 2 ps nstvout = 1000 ; save velocities every 2 ps nstxtcout = 1000 ; xtc compressed trajectory output every 2 ps nstenergy = 1000 ; save energies every 2 ps nstlog = 1000 ; update log file every 2 ps ; Bond parameters continuation= yes ; Restarting after NPT constraint_algorithm = lincs; holonomic constraints constraints = all-bonds ; all bonds (even heavy atom-H bonds) constrained lincs_iter = 1 ; accuracy of LINCS lincs_order = 4 ; also related to accuracy ; Neighborsearching ns_type = grid ; search neighboring grid cells nstlist = 5 ; 10 fs rlist = 1.0 ; short-range neighborlist cutoff (in nm) rcoulomb= 1.0 ; short-range electrostatic cutoff (in nm) rvdw= 1.0 ; short-range van der Waals cutoff (in nm) ; Electrostatics coulombtype = PME ; Particle Mesh Ewald for long-range electrostatics pme_order = 4 ; cubic interpolation fourierspacing = 0.16 ; grid spacing for FFT ; Temperature coupling is on tcoupl = V-rescale ; modified Berendsen thermostat tc-grps = Protein Non-Protein ; two coupling groups - more accurate tau_t = 0.1 0.1 ; time constant, in ps ref_t = 320 320 ; reference temperature, one for each group, in K ; Pressure coupling is on pcoupl = Parrinello-Rahman ; Pressure coupling on in NPT pcoupltype = isotropic ; uniform scaling of box vectors tau_p = 2.0 ; time constant, in ps ref_p = 1.0 ; reference pressure, in bar compressibility = 4.5e-5; isothermal compressibility of water, bar^-1 ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr= EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = no; Velocity generation is off -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] trjconv -pbc nojump across multiple trajectories
Hi all, I am planning to run a 100ns simulation by continuing a simulation in increments of 1ns. After each round, analyses are performed and the trajectory scrapped. One of the analysis I need to do is mean square displacement; for this I need a continuous trajectory as provided by trjconv -pbc nojump. The issue that I have is how to make sure that the trajectory stays continuous across the different time increments. By reading the manual of trjconv, I found out that the starting configuration for this procedure is taken from the structure file, if one is supplied, otherwise it is the first frame. I thought of using the .cpt file from the previous step as structure file (trjconv -s) for trjconv, but .cpt is not one of the formats accepted for that option (only tpr, tpb, tpa, gro, g96, pdb). The other issue is that as the dumping of snapshots is performed frequently (every 100 timesteps) I would be using double precision, so I would like to use a double-precision structure to do the nojump step as well. A cpt file seems to be perfectly suited for this, as well as a single-frame trr file, however neither of them is accepted as a structure file for trjconv. My questions are 2: 1. what would be the best way to ensure a continuous trajectory across multiple fragments of a trajectory, using double precision? 2. is there a specific reason why a cpt is not accepted as a structure file for trjconv -s? Thanks for your comments! Regards, Pablo -- Dr. Pablo Englebienne Postdoctoral Researcher *TU Delft / 3mE / Process Energy* /Engineering Thermodynamics (ETh) group/ Building 46 Leeghwaterstraat 44, room 030 2628 CA Delft The Netherlands *T* +31 (0)15 27 86662 tel:+31152786662 *E* p.englebie...@tudelft.nl mailto:p.englebie...@tudelft.nl -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Re: trjconv -pbc nojump across multiple trajectories
Thanks for your reply, Tsjerk. Indeed, trjcat is a good option. I could first run trjconv -pbc nojump on a trajectory, then keep the last frame as a trr and stitch it with the following trajectory; trjconv then takes the first frame as reference for the nojump. Regarding this, what would be the easiest way to get the last frame out of the trajectory? I can think of using the output of gmxcheck to get the time of the last frame and then use it as the argument to trjconv -b, but is there a simpler way? Regards, Pablo Dr. Pablo Englebienne Postdoctoral Researcher *TU Delft / 3mE / Process Energy* /Engineering Thermodynamics (ETh) group/ Building 46 Leeghwaterstraat 44, room 030 2628 CA Delft The Netherlands *T* +31 (0)15 27 86662 tel:+31152786662 *E* p.englebie...@tudelft.nl mailto:p.englebie...@tudelft.nl On 11/22/2012 04:20 PM, gmx-users-requ...@gromacs.org wrote: Send gmx-users mailing list submissions to gmx-users@gromacs.org To subscribe or unsubscribe via the World Wide Web, visit http://lists.gromacs.org/mailman/listinfo/gmx-users or, via email, send a message with subject or body 'help' to gmx-users-requ...@gromacs.org You can reach the person managing the list at gmx-users-ow...@gromacs.org When replying, please edit your Subject line so it is more specific than Re: Contents of gmx-users digest... Today's Topics: 1. Re: Re,i don't know how can i determine emtol (Ivan Gladich) 2. top2psf with Amber99SB-ILDN (Steven Neumann) 3. Re: top2psf with Amber99SB-ILDN (francesco oteri) 4. Re: top2psf with Amber99SB-ILDN (Steven Neumann) 5. Bonded parametrs for CG (Steven Neumann) 6. trjconv -pbc nojump across multiple trajectories (Pablo Englebienne) 7. Re: trjconv -pbc nojump across multiple trajectories (Tsjerk Wassenaar) 8. On the usage of SD integrator as the thermostat (James Starlight) -- Message: 1 Date: Thu, 22 Nov 2012 12:17:50 +0100 From: Ivan Gladich ivan.glad...@marge.uochb.cas.cz Subject: Re: [gmx-users] Re,i don't know how can i determine emtol To: Discussion list for GROMACS users gmx-users@gromacs.org Cc: Ali Alizadeh ali.alizadehmoja...@gmail.com Message-ID: 50ae09de.5040...@marge.uochb.cas.cz Content-Type: text/plain; charset=UTF-8; format=flowed On 11/21/2012 06:02 PM, Ali Alizadeh wrote: 1- In your opinion, Can i simulate that system? In my (humble) opinion: 1)Of course you can simulate that system...however I doubt that, without starting from the exact initial configuration with the exactly same set-up, you can get the same results (i.e. see the nucleation). The onset of ice nucleation is a random process and requires very long simulation (the paper that you posted was analysing micro-second trajectories!!). There is the risk that you could try several different initial configurations at several temperature without getting anything. However, read carefully that paper, I do not remember all the details. If you are interested in ice crystal growth, I would suggest to start with an initial water/ice system: at temperature below the melting one you will see formation of new ice starting from the initial ice matrix. I do not know how can ice crystal in initial step, i study papers and i know ice , crystals formation and dissociation of crystals is very difficult, and we can use this method(as you said above) for melting or freezing(if i say correct) My problem is this case, how can i have a crystal in initial step, You have three options 1) Surfing the web searching for some structure. 2)Try to implement the proton disordering algorithm of Buch et al., (V. Buch, P. Sandler and J. Sadlej, J. Phys. Chem. B, 1998, 102, 8641âEUR8653) which specifies orientations of water molecules such that ice Ih BernalâEURFowler constraints for each molecule are satisfied. (J. D. Bernal and R. H. Fowler, J. Chem. Phys., 1933, 1, 515âEUR548) 3) Someone kindly offer you his own structure, e.g. http://marge.uochb.cas.cz/~gladich/Teaching.html (PS: This is the structure at O K...you have to heat it to the desired temperature) Good Luck ivan There are several works in literature on this. 2- How can i use rigid TIP4P model of water? 2) I would rather use other water model that have been explicitly tested for Ice (e.g. TIP4P/2005, TIP4P/Ice, TIP5P-Ew, NE6) I think that i solved it, thank you, Best Ivan Sincerely Ali Alizadeh On 11/21/2012 09:43 AM, Ivan Gladich wrote: Dear Ali the paper that you are citing is using a rigid TIP4P water model As far as I know, emtol is relevant only for minimization or molecular dynamics with shell particle or flexible constraints. Therefore, as Justin told you, the emtol value should be irrelevant. Concerning this paper, I would like to warn you that h_omogeneous ice nucleation from bulk water_ with explicit water molecule is very rare event... It's depends from
[gmx-users] g(r) does not go to 1 at long r -- bug in g_rdf?
Hi, I tried to calculate the radial distribution functions for a simple system: a 5nm a side cubic box with 10 Ne atoms and 10 Ar atoms, simulated for 100ns in NVT @ 300K. I was expecting to get an RDF with a peak, stabilizing to 1.0 at long distances. This was the case for the Ne-Ar RDF, but not for the Ne-Ne or Ar-Ar RDFs, which stabilize to about 0.9. I believe this is due to a problem in the normalization of the histograms with respect to the number of pairs available: there are N*N pairs for the Ne-Ar, while N*(N-1) for the Ne-Ne and Ar-Ar case. Did somebody else find an issue like this? I think that the issue may become not evident for a relatively large system, as N*N ~ N*(N-1) for large N. I put the relevant files if somebody wishes to reproduce it here: https://gist.github.com/4085292 I'll appreciate input on this and I can also file a bug if deemed necessary. Take care, Pablo -- Dr. Pablo Englebienne Postdoctoral Researcher *TU Delft / 3mE / Process Energy* /Engineering Thermodynamics (ETh) group/ Building 46 Leeghwaterstraat 44, room 030 2628 CA Delft The Netherlands *T* +31 (0)15 27 86662 tel:+31152786662 *E* p.englebie...@tudelft.nl mailto:p.englebie...@tudelft.nl -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! * Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Re: Difficulty building a topology for a synthetic, branched PEG-peptide molecule [SOLVED]
That's a long bond. Does your reference length in specbond.dat suit it? IIRC there should be some evidence in the output of the special bond being formed if it actually is. If not, your symptoms are probably related Hi Mark, indeed, I think that was part of the problem. pdb2gmx indeed outputs a message when a specbond.dat rule is matched and a bond formed. In case it helps to someone else or for reference purposes, I finally managed to solve the issue. Some of the lessons I learned: - make sure that the residues/atoms in specbond.dat were correct (I had defined a number of residues for termini and mid-chain PEG and connectors, and I got them confused at some point, so not all of them were being recognized properly). The dangling bond at at least one of the terminal ends given by pdb2gmx is most likely due to this and/or the protonation state of the residues connecting the fragments - each fragment is assigned a different chain letter in the source PDB file; in my case, for a system like [N-(peptide1)-C]-[N-(PEG)-C]-[N-(peptide2)-Lys-C] | NZ | [C-(PEG)-N]-[C-(peptide3)-N] each fragment in square brackets is assigned a different chain name in [A-E] in the PDB file: ATOM123 N ARG A 1 74.024 13.299 50.237 1.00 0.00 N1+ ^ - always list the atoms within a chain in an N-to-C direction. This means that the main branch is defined first, and then the branching point is defined in an N-to-C direction, even if it is counterintuitive by the way they are connected. specbond.dat takes good care of setting up the connection in all cases (as long as they are well defined...). - the PEG residues need to be defined as type Other in residues.dat - call pdb2gmx to manually assign the termini and Lys protonation states manually, and to merge the chains into a single molecule: pdb2gmx-ter -f substrate.pdb-chainsep interactive-lys - the internal termini (i.e., the ones that are peptide termini but are connected to a PEG chain) need to be given a protonation state of None, while the real termini can be assigned as charged or neutral, depending on the conditions Oh, and well done for constructing a good question. You would likely not have gotten anywhere giving less detail :) Thanks, it actually helped putting everything in writing, as it pointed out the few things that I hadn't yet looked at in detail... Pablo Englebienne, PhD Dept. of Biomedical Engineering Dept. of Chemistry and Chemical Engineering Institute for Complex Molecular Systems (ICMS) Eindhoven University of Technology, TU/e PO Box 513, HG -1.26 5600 MB Eindhoven, The Netherlands Tel +31 40 247 5349 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Difficulty building a topology for a synthetic branched PEG-peptide molecule
Hi all, I am trying to build a topology for a synthetic molecule that consists of peptides connected by oligoethyleneglycol (I'll call it PEG) linkers terminated with an amine and a carboxylic acid: -NH2-CH2-[CH2-O-CH2]n-CH2-C(=O)- The system looks like this: N-(PEG)-C-N-(peptide2)-Lys(C-ter)-NZ-C-(PEG)-N So: * the C-terminus of a PEG linker is attached to the N terminus of the peptide * the terminal Lys on the peptide is attached to the C-terminus of a PEG linker I was able to successfully build a topology for this molecule by: * defining appropriate residues (for the PEG chains and the Lys with a PEG attached on the NZ) in a local copy of the forcefield file, adding the residues' topologies to aminoacids.rtp * using the specbond.dat file to define the bond between the NZ in the Lys and the PEG linker * adding the residues to residuetypes.dat with a Protein type * calling pdb2gmx with the -ter option to assign the protonation states Now, I need to extend the topology to a molecule like this one: N-(peptide1)-C-N-(PEG)-C-N-(peptide2)-Lys(C-ter)-NZ-C-(PEG)-N-C-(peptide3)-N The difficulty with this molecule is that it has 2 N-termini and a single C-terminus (in the Lys with the PEG attached to the NZ sidechain). pdb2gmx recognizes the whole molecule as a peptide, but treats the last residue as a C-terminus, when it actually is an N-terminus. I found in the description for specbond.dat (http://www.gromacs.org/Documentation/File_Formats/specbond.dat) that for a branched peptide the -chainsep option of pdb2gmx can be used, so I started to work on that. I split the molecules in 2 chains like this: N-(peptide1)-C-N-(PEG)-C-N-(peptide2)-Lys(C-ter)-NZ- || -C-(PEG)-N-C-(peptide3)-N I reversed the order of the residues in the second chain, so that the residues are in N-to-C order. With this, pdb2gmx recognizes the proper termini when called as: $ pdb2gmx -f substrate_edited-reversed.pdb -ter -chainsep interactive I tried setting the protonation of the termini as charged for the real termini and None for the artificial one (the one that should be handled by specbond.dat), but I get an error message: ---[pdb2gmx output]--- [...] Splitting PDB chains based on TER records or changing chain id. Merge chain ending with residue LYSS27 (chain id ' ', atom 71 NZ) with chain starting with residue GLU28 (chain id 'p', atom 308 OE2)? [n/y] n Merged 1 chains into one molecule definition There are 2 chains and 0 blocks of water and 47 residues with 303 atoms chain #res #atoms 1 ' '27198 2 'p'17105 [...] Identified residue ARG1 as a starting terminus. Identified residue LYSS27 as a ending terminus. 9 out of 9 lines of specbond.dat converted successfully Special Atom Distance matrix: PEA118 CG1125 LYSS27 NZ198 2.762 Select start terminus type for ARG-1 0: NH3+ 1: NH2 2: None 0 Start terminus ARG-1: NH3+ Select end terminus type for LYSS-27 0: COO- 1: COOH 2: None 0 End terminus LYSS-27: COO- [...] Identified residue GLU28 as a starting terminus. Identified residue PEA47 as a ending terminus. 9 out of 9 lines of specbond.dat converted successfully Select start terminus type for GLU-28 0: NH3+ 1: NH2 2: None 0 Start terminus GLU-28: NH3+ Select end terminus type for PEA-47 0: COO- 1: COOH 2: None 2 End terminus PEA-47: None [...] --- Program pdb2gmx, VERSION 4.5.4 Source code file: pdb2top.c, line: 1035 Fatal error: There is a dangling bond at at least one of the terminal ends. Select a proper terminal entry. For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors --- ---[pdb2gmx output]--- I think that the problem might stem from the fact that the C-terminus in the second chain is not a real peptide; I changed residuetypes.dat to have the PEG residues as Other, which causes pdb2gmx to recognize the last aminoacid as a C-terminus, but treating it as None yields the same error. Unfortunately, there is nothing about this error in http://www.gromacs.org/Documentation/Errors . Any suggestions on how to make this work will be greatly appreciated! -- Pablo Englebienne, PhD Dept. of Biomedical Engineering Dept. of Chemistry and Chemical Engineering Institute for Complex Molecular Systems (ICMS) Eindhoven University of Technology, TU/e PO Box 513, HG -1.26 5600 MB Eindhoven, The Netherlands Tel +31 40 247 5349 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Error compiling Gromacs 4.5.4: relocation R_X86_64_32 against `a local symbol' can not be used when making a shared object; recompile with -fPIC
Hi all, I'm trying to compile release 4.5.4 on a system that has been running every release since 4.0.4 without a problem. Even 4.5.3 compiled fine with the following configure: LDFLAGS=-L/cvos/shared/apps/fftw/gcc/64/3.2/lib CPPFLAGS=-I/cvos/shared/apps/fftw/gcc/64/3.2/include ./configure --prefix=$HOME/software The LDFLAGS and CPPFLAGS specify the (non-standard) location of the FFTW libraries and headers. Configure succeeds in creating the Makefiles, but when running make it aborts at this point: cc -shared .libs/calcmu.o .libs/calcvir.o .libs/constr.o .libs/coupling.o .libs/domdec.o .libs/domdec_box.o .libs/domdec_con.o .libs/domdec_network.o .libs/domdec_setup.o .libs/domdec_top.o .libs/ebin.o .libs/edsam.o .libs/ewald.o .libs/force.o .libs/forcerec.o .libs/ghat.o .libs/init.o .libs/mdatom.o .libs/mdebin.o .libs/minimize.o .libs/mvxvf.o .libs/ns.o .libs/nsgrid.o .libs/perf_est.o .libs/genborn.o .libs/genborn_sse2_single.o .libs/genborn_sse2_double.o .libs/genborn_allvsall.o .libs/genborn_allvsall_sse2_single.o .libs/genborn_allvsall_sse2_double.o .libs/gmx_qhop_parm.o .libs/gmx_qhop_xml.o .libs/groupcoord.o .libs/pme.o .libs/pme_pp.o .libs/pppm.o .libs/partdec.o .libs/pull.o .libs/pullutil.o .libs/rf_util.o .libs/shakef.o .libs/sim_util.o .libs/shellfc.o .libs/stat.o .libs/tables.o .libs/tgroup.o .libs/tpi.o .libs/update.o .libs/vcm.o .libs/vsite.o .libs/wall.o .libs/wnblist.o .libs/csettle.o .libs/clincs.o .libs/qmmm.o .libs/gmx_fft.o .libs/gmx_parallel_3dfft.o .libs/fft5d.o .libs/gmx_wallcycle.o .libs/qm_gaussian.o .libs/qm_mopac.o .libs/qm_gamess.o .libs/gmx_fft_fftw2.o .libs/gmx_fft_fftw3.o .libs/gmx_fft_fftpack.o .libs/gmx_fft_mkl.o .libs/qm_orca.o .libs/mdebin_bar.o -Wl,--rpath -Wl,/home/penglebie/downloads/gromacs-4.5.4/src/gmxlib/.libs -Wl,--rpath -Wl,/home/penglebie/software/lib /cvos/shared/apps/fftw/gcc/64/3.2/lib/libfftw3f.a -lxml2 -L/cvos/shared/apps/fftw/gcc/64/3.2/lib ../gmxlib/.libs/libgmx.so -lnsl -lm -msse2 -pthread -Wl,-soname -Wl,libmd.so.6 -o .libs/libmd.so.6.0.0 /usr/bin/ld: /cvos/shared/apps/fftw/gcc/64/3.2/lib/libfftw3f.a(plan-many-dft-r2c.o): relocation R_X86_64_32 against `a local symbol' can not be used when making a shared object; recompile with -fPIC /cvos/shared/apps/fftw/gcc/64/3.2/lib/libfftw3f.a: could not read symbols: Bad value collect2: ld returned 1 exit status make[3]: *** [libmd.la] Error 1 make[3]: Leaving directory `/home/penglebie/downloads/gromacs-4.5.4/src/mdlib' make[2]: *** [all-recursive] Error 1 make[2]: Leaving directory `/home/penglebie/downloads/gromacs-4.5.4/src' make[1]: *** [all] Error 2 make[1]: Leaving directory `/home/penglebie/downloads/gromacs-4.5.4/src' make: *** [all-recursive] Error 1 I see that recently (http://lists.gromacs.org/pipermail/gmx-users/2011-April/059919.html) another user encountered the same problem but this time with version 4.5.3; in my case 4.5.3 compiles fine, the only issue is with 4.5.4. The system is running Scientific Linux 5.5. $ uname -a Linux ST-HPC-Main 2.6.18-128.7.1.el5 #1 SMP Mon Aug 24 08:12:52 EDT 2009 x86_64 x86_64 x86_64 GNU/Linux I am puzzled as to why it doesn't work in 4.5.4 but did until the previous release. Did something change in this respect? Regards, Pablo -- Pablo Englebienne, PhD Dept. of Biomedical Engineering Dept. of Chemistry and Chemical Engineering Institute for Complex Molecular Systems (ICMS) Eindhoven University of Technology, TU/e PO Box 513, HG -1.26 5600 MB Eindhoven, The Netherlands Tel +31 40 247 5349 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/Search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
[gmx-users] Continuous trajectory from trjconv?
Hi, I made an NpT MD simulation of 10 copies of the same small molecule in CHCl3 (dodecahedron unit cell) for 20 ns. I want to visualize the evolution of the system and the distance between the molecules. I tried to visualize the system, but I can't get a continuous trajectory (i.e., without jumps) without the solutes diffusing out of the simulation box. The same thing happens for the inter-molecule distance, it oscillates wildly between a reasonable value (say, 0.5 nm) and 1/2 the size of the periodic cell. When I use trjconv with the -pbc nojump -ur compact options, the molecules diffuse out of the box, in a continuous trajectory, until they are way further apart than the size of the cell. When visualizing it in VMD, however, it turns out that not all molecules are isolated, but instead some are interacting with the periodic image of another molecule. I used pbc wrap from the PBCTools module within VMD and that brings all the molecules into a single cell, but then the trajectory is not continuous, with compounds jumping around the edges of the unit cell. It is then not straightforward to see the interaction among different molecules, as in some cases they are in opposite edges of the cell. If I use the -pbc mol -ur compact option, all molecules stay within a single unit cell, but there are jumps across the border of the unit cell as in the above case. I tried also using the -center, -fit rot+trans and -fit progressive options (with groups containing either all molecules or a single residue) but this ultimately gave the same results. Would there be something else I could try to visualize the simulation so that the central unit cell would contain all residues in a continuous way? Thanks in advance for reading, and more for replying! Take care, Pablo -- Pablo Englebienne, PhD Institute of Complex Molecular Systems (ICMS) Eindhoven University of Technology, TU/e PO Box 513, HG -1.26 5600 MB Eindhoven, The Netherlands Tel +31 40 247 5349 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] Re: Continuous trajectory from trjconv?
Thanks Tsjerk, I figured as much, but wanted to make sure I wasn't overlooking something... Take care, Pablo Hi Pablo, You want to mutually exclusive things. That is by definition impossible. Sorry, Tsjerk On Mon, Mar 1, 2010 at 11:02 AM, Pablo Englebienne p.englebie...@tue.nl wrote: Hi, I made an NpT MD simulation of 10 copies of the same small molecule in CHCl3 (dodecahedron unit cell) for 20 ns. I want to visualize the evolution of the system and the distance between the molecules. I tried to visualize the system, but I can't get a continuous trajectory (i.e., without jumps) without the solutes diffusing out of the simulation box. The same thing happens for the inter-molecule distance, it oscillates wildly between a reasonable value (say, 0.5 nm) and 1/2 the size of the periodic cell. When I use trjconv with the -pbc nojump -ur compact options, the molecules diffuse out of the box, in a continuous trajectory, until they are way further apart than the size of the cell. When visualizing it in VMD, however, it turns out that not all molecules are isolated, but instead some are interacting with the periodic image of another molecule. I used pbc wrap from the PBCTools module within VMD and that brings all the molecules into a single cell, but then the trajectory is not continuous, with compounds jumping around the edges of the unit cell. It is then not straightforward to see the interaction among different molecules, as in some cases they are in opposite edges of the cell. If I use the -pbc mol -ur compact option, all molecules stay within a single unit cell, but there are jumps across the border of the unit cell as in the above case. I tried also using the -center, -fit rot+trans and -fit progressive options (with groups containing either all molecules or a single residue) but this ultimately gave the same results. Would there be something else I could try to visualize the simulation so that the central unit cell would contain all residues in a continuous way? Pablo Englebienne, PhD Institute of Complex Molecular Systems (ICMS) Eindhoven University of Technology, TU/e PO Box 513, HG -1.26 5600 MB Eindhoven, The Netherlands Tel +31 40 247 5349 -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] How to tune number of CPUs for a run?
Hi all,
I'm having some trouble running simulations with increasing number of
CPUs. What parameters should I modify to make sure that the simulation
would run with a specific number of processors? Or, having access to a
large number of processors, how to select the number of CPUs to request?
Besides this, should the PP/PME reported by grompp always fall in the
range 0.25-0.33? What if it is lower (e.g., 0.16)?
I'm attaching an mdrun logfile of a failed run.
Thanks for suggestions,
Pablo
--
Pablo Englebienne, PhD
Institute of Complex Molecular Systems (ICMS)
Eindhoven University of Technology, TU/e
PO Box 513, HG -1.26
5600 MB Eindhoven, The Netherlands
Tel +31 40 247 5349
Log file opened on Mon Nov 2 18:23:16 2009
Host: node052 pid: 22760 nodeid: 0 nnodes: 16
The Gromacs distribution was built Thu Oct 29 14:19:59 CET 2009 by
pengle...@st-hpc-main (Linux 2.6.18-128.7.1.el5 x86_64)
:-) G R O M A C S (-:
Good gRace! Old Maple Actually Chews Slate
:-) VERSION 4.0.5 (-:
Written by David van der Spoel, Erik Lindahl, Berk Hess, and others.
Copyright (c) 1991-2000, University of Groningen, The Netherlands.
Copyright (c) 2001-2008, The GROMACS development team,
check out http://www.gromacs.org for more information.
This program is free software; you can redistribute it and/or
modify it under the terms of the GNU General Public License
as published by the Free Software Foundation; either version 2
of the License, or (at your option) any later version.
:-) /home/penglebie/software/bin/mdrun_openmpi (double precision) (-:
PLEASE READ AND CITE THE FOLLOWING REFERENCE
B. Hess and C. Kutzner and D. van der Spoel and E. Lindahl
GROMACS 4: Algorithms for highly efficient, load-balanced, and scalable
molecular simulation
J. Chem. Theory Comput. 4 (2008) pp. 435-447
--- Thank You ---
PLEASE READ AND CITE THE FOLLOWING REFERENCE
D. van der Spoel, E. Lindahl, B. Hess, G. Groenhof, A. E. Mark and H. J. C.
Berendsen
GROMACS: Fast, Flexible and Free
J. Comp. Chem. 26 (2005) pp. 1701-1719
--- Thank You ---
PLEASE READ AND CITE THE FOLLOWING REFERENCE
E. Lindahl and B. Hess and D. van der Spoel
GROMACS 3.0: A package for molecular simulation and trajectory analysis
J. Mol. Mod. 7 (2001) pp. 306-317
--- Thank You ---
PLEASE READ AND CITE THE FOLLOWING REFERENCE
H. J. C. Berendsen, D. van der Spoel and R. van Drunen
GROMACS: A message-passing parallel molecular dynamics implementation
Comp. Phys. Comm. 91 (1995) pp. 43-56
--- Thank You ---
parameters of the run:
integrator = md
nsteps = 5
init_step= 0
ns_type = Grid
nstlist = 5
ndelta = 2
nstcomm = 1
comm_mode= Linear
nstlog = 1000
nstxout = 1000
nstvout = 1000
nstfout = 0
nstenergy= 1000
nstxtcout= 0
init_t = 0
delta_t = 0.002
xtcprec = 1000
nkx = 40
nky = 40
nkz = 40
pme_order= 4
ewald_rtol = 1e-05
ewald_geometry = 0
epsilon_surface = 0
optimize_fft = FALSE
ePBC = xyz
bPeriodicMols= FALSE
bContinuation= TRUE
bShakeSOR= FALSE
etc = V-rescale
epc = Parrinello-Rahman
epctype = Isotropic
tau_p= 5
ref_p (3x3):
ref_p[0]={ 1.0e+00, 0.0e+00, 0.0e+00}
ref_p[1]={ 0.0e+00, 1.0e+00, 0.0e+00}
ref_p[2]={ 0.0e+00, 0.0e+00, 1.0e+00}
compress (3x3):
compress[0]={ 1.0e-04, 0.0e+00, 0.0e+00}
compress[1]={ 0.0e+00, 1.0e-04, 0.0e+00}
compress[2]={ 0.0e+00, 0.0e+00, 1.0e-04}
refcoord_scaling = No
posres_com (3):
posres_com[0]= 0.0e+00
posres_com[1]= 0.0e+00
posres_com[2]= 0.0e+00
posres_comB (3):
posres_comB[0]= 0.0e+00
posres_comB[1]= 0.0e+00
posres_comB[2]= 0.0e+00
andersen_seed= 815131
rlist= 1.4
rtpi = 0.05
coulombtype = PME
rcoulomb_switch = 0
rcoulomb = 1.4
vdwtype = Cut-off
rvdw_switch = 0
rvdw = 1.4
epsilon_r= 1
epsilon_rf = 1
tabext = 1
implicit_solvent = No
gb_algorithm = Still
[gmx-users] Re: Chloroform (CHCl3) solvent box for G53a5 force field
in potential energy (-4000 kJ; RMSD 1000) - tau_p = 1.0 == increase of Density to 1850 in the first 10ps, then stable with fluctuations of ~34 in the remaining 90ps; potential energy with huge fluctuations - tau_p = 0.2 (as used in Mol. Phys. 1994, 83, 381) == LINCS warnings and posterior crash: Step 23, time 0.046 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.079726, max 0.094860 (between atoms 672 and 673) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length 242243 36.20.2504 0.1919 0.1758 [...] --- Program mdrun_mpi_d, VERSION 4.0.5 Source code file: nsgrid.c, line: 357 Range checking error: Explanation: During neighborsearching, we assign each particle to a grid based on its coordinates. If your system contains collisions or parameter errors that give particles very high velocities you might end up with some coordinates being +-Infinity or NaN (not-a-number). Obviously, we cannot put these on a grid, so this is usually where we detect those errors. Make sure your system is properly energy-minimized and that the potential energy seems reasonable before trying again. Variable ci has value -2147483631. It should have been within [ 0 .. 84 ] --- I'll appreciate pointers as to what to try next! Regards, Pablo -- Pablo Englebienne, PhD Institute of Complex Molecular Systems (ICMS) Eindhoven University of Technology, TU/e PO Box 513, HG -1.26 5600 MB Eindhoven, The Netherlands Tel +31 40 247 5349 ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] Re: Chloroform (CHCl3) solvent box for G53a5 force field
Following up on the previous message, I noticed that the topology I previously sent (including 10 bonds) works for a minimization, but not for an MD simulation. grompp issues the following warning: WARNING 1 [file topol.top, line 29]: Molecule type 'CHCL3' has 10 constraints. For stability and efficiency there should not be more constraints than internal number of degrees of freedom: 9. I therefore used the following .itp file, with only the C-H and C-Cl bonds: --[chcl3.itp]-- [ moleculetype ] ; Namenrexcl CHCL3 1 [ atoms ] ; nr type resnr residue atom cgnr charge mass typeBchargeB massB 1 HCHL 1 CHCL3 HChL 1 0.082 1.008 2 CCHL 1 CHCL3 CChL 1 0.179 12.011 3 CLCHL 1 CHCL3 CLCh1 1 -0.087 35.453 4 CLCHL 1 CHCL3 CLCh2 1 -0.087 35.453 5 CLCHL 1 CHCL3 CLCh3 1 -0.087 35.453 [ bonds ] ; aiaj functc0c1c2c3 1 2 2gb_39 2 3 2gb_40 2 4 2gb_40 2 5 2gb_40 [ angles ] ; aiajak functc0c1 c2c3 1 2 3 2 ga_43 1 2 4 2 ga_43 1 2 5 2 ga_43 3 2 4 2 ga_44 3 2 5 2 ga_44 4 2 5 2 ga_44 [ exclusions ] 1 2 3 4 5 2 1 3 4 5 3 1 2 4 5 4 1 2 3 5 5 1 2 3 4 --[chcl3.itp]-- These are the steps I took: - build a box with 216 CHCL3 molecules with genbox - adjusted the density to 1479 with editconf - minimized the box to F100 (steep, 5 steps) - equilibrated NVT (position restrained C), 100ps, 300K, tau_t 0.1 - equilibrated NPT (position restrained C), 100ps, 300K, 1bar, tau_t 0.1, tau_p 2.0, compressibility 1e-4 (from CRC handbook) Until here, everything looks decent, except for relatively large fluctuations in T (RMSD ~9K) and density (RMSD ~10kg m-3). I then performed an unconstrained MD, 1ns, otherwise identical parameters to NPT equilibration. Temp (300K) and density (1450 kg m-3) stable, but fluctuating (RMSD 8 and 20 respectively). In order to compare the results with Tironi and van Gunsteren, Molecular Physics 1994, 83, 381, who used the same GROMOS parameters: Epot = -28.6 +/- 0.3 kJ/mol, density 1520 +/- 12 kg m-3. This is the output I get from g_energy: Energy Average RMSD Fluct. Drift Tot-Drift --- Potential -4720.16119.141118.978 -0.0215895 -21.5896 Therefore molar Epot = -(-4720)/216 = 21.9 kJ/mol. What factors could be accountable for the decrease wrt the reported value of 28.4? Experimental deltaHv is 31.4 kJ/mol (same reference). Would there be any other parameters I should check before using this solvent box in production runs? Regards, Pablo -- Pablo Englebienne, PhD Institute of Complex Molecular Systems (ICMS) Eindhoven University of Technology, TU/e PO Box 513, HG -1.26 5600 MB Eindhoven, The Netherlands Tel +31 40 247 5349 ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] Re: Chloroform (CHCl3) solvent box for G53a5 force field
OK, I started over with the CHCl3 box from scratch. I prepared the following itp file from the CHCL3 parameters in ffG53a5.rtp: ---[chcl3.itp]--- [ moleculetype ] ; Namenrexcl CHCL3 3 [ atoms ] ; nr type resnr residue atom cgnr charge mass 1 HCHL 1 CHCL3 HChL 1 0.082 1.008 2 CCHL 1 CHCL3 CChL 1 0.179 12.011 3 CLCHL 1 CHCL3 CLCh1 1 -0.087 35.453 4 CLCHL 1 CHCL3 CLCh2 1 -0.087 35.453 5 CLCHL 1 CHCL3 CLCh3 1 -0.087 35.453 [ bonds ] ; aiaj funct 1 2 2gb_39 2 3 2gb_40 2 4 2gb_40 2 5 2gb_40 1 3 2gb_47 1 4 2gb_47 1 5 2gb_47 3 4 2gb_48 3 5 2gb_48 4 5 2gb_48 [ angles ] ; aiajak funct 1 2 3 2 ga_43 1 2 4 2 ga_43 1 2 5 2 ga_43 3 2 4 2 ga_44 3 2 5 2 ga_44 4 2 5 2 ga_44 ---[chcl3.itp]--- I noticed that G53a5 includes 4 types of bond stretching terms specific for CHCl3 (C-Cl, C-H, H-Cl and Cl-Cl), therefore I specified all of them. Should all of these terms be harmonic bonds (function type 2) or some (e.g., the H-Cl and Cl-Cl terms) should be type 6 (as described in section 5.4 of the manual)? I tried with both types and I get the same minimized structure with the following topology and mdp files: ---[topol.top]--- ; Include forcefield parameters #include ffG53a5.itp #include chcl3.itp [ system ] ; Name Chloroform [ molecules ] ; Compound#mols CHCL3 1 ---[topol.top]--- ---[em.mdp]--- integrator = cg nsteps = 5 ; ; Energy minimizing stuff ; emtol = 1e-5 emstep = 0.01 ; ; Electrostatics ; coulombtype = cut-off rcoulomb= 0 ns_type = simple nstlist = 0 rlist = 0 ; ; vdW ; rvdw= 0 ; ; PBC ; pbc = no ---[em.mdp]--- In the minimized structure, not all C-Cl distances are equivalent, although the minimization converges OK. Any comments? ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] Chloroform (CHCl3) solvent box for G53a5 force field
Hi all, I'm trying to simulate a small molecule in a chloroform box using the GROMOS G53a5 forcefield. I realized that the parameters for the solvent are present in the ffG53a5.rtp file, however I could not find a CHCl3 solvent box included in GROMACS. I did find, however, a CHCl3 solvent box equilibrated by PeiQuan Chen (http://lists.gromacs.org/pipermail/gmx-users/2003-May/005572.html) at http://www.gromacs.org/index.php?title=Download_%26_Installation/User_contributions/Molecule_topologies . I also saw mention of a box for the Amber forcefield, that is now included in AmberTools (amber10/dat/solvents/cform/cform.pdb and chcl3_equil.pdb.1). I took this one and I'm equilibrating it to use it later. In the mean time, I wanted to know if I overlooked something, and there is a CHCl3 box to use with the GROMOS forcefield? Thanks! -- Pablo Englebienne, PhD Institute of Complex Molecular Systems (ICMS) Eindhoven University of Technology, TU/e PO Box 513, HG -1.26 5600 MB Eindhoven, The Netherlands Tel +31 40 247 5349 ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] Re: Chloroform (CHCl3) solvent box for G53a5 force field
I'm trying to simulate a small molecule in a chloroform box using the GROMOS G53a5 forcefield. I realized that the parameters for the solvent are present in the ffG53a5.rtp file, however I could not find a CHCl3 solvent box included in GROMACS. I did find, however, a CHCl3 solvent box equilibrated by PeiQuan Chen (http://lists.gromacs.org/pipermail/gmx-users/2003-May/005572.html) at http://www.gromacs.org/index.php?title=Download_%26_Installation/User_contributions/Molecule_topologies . I also saw mention of a box for the Amber forcefield, that is now included in AmberTools (amber10/dat/solvents/cform/cform.pdb and chcl3_equil.pdb.1). I took this one and I'm equilibrating it to use it later. In the mean time, I wanted to know if I overlooked something, and there is a CHCl3 box to use with the GROMOS forcefield? Not one that is officially distributed. If it was, it would be the /share/gromacs/top subdirectory with other solvent topologies and structures. Probably your best bet is to use the one in the User Contributions section, unless you feel the need to create your own and start from scratch. Thanks for confirming this, Justin. I decided not to use the user-contributed CHCl3 box because the topology is not consistent with the GROMOS atom types: the CH is united atom (although the mass is 12.01100, it should probably be 13.01900?), while in G53a5 there are parameters for C, H and Cl. Has someone ever used this box successfully? I was able to equilibrate the box in NVT (100-200ps gives a stable simulation), although it exploded at constant pressure. After looking at some references on the mailing lists (both GMX and AMBER), I tried increasing tau_p from 2.0 to 5.0 and that yielded a stable 100 ps simulation, although the system later (continuing for further ~150 ps) started to oscillate wildly in temperature and pressure. What is the effect of increasing tau_p? Besides making the dynamics stable, would changing its value affect the outcome of the simulation in any other way? Thanks again! ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
[gmx-users] Improper dihedrals on planar rings
Hi all, I have a question about the setting of improper dihedrals. I have a pyrimidinone I want to use with the G53a5 force field, but I can't get the ring to stay flat. First question is: do all atoms in the planar ring need to have 3 subsituents? I.e., should an explicit atom of type HC be present on the CHs within the ring, or could they be united atoms (e.g., non-bonded type CR1)? Second question: how many parameters are required to keep the ring planar? I initially though that each bond within the ring should have an improper torsion defined, with the bonds in the J-K atoms in the torsion. In addition, every center bound to 3 other atoms (e.g., C=O, N-H) needs an additional term, with the central atom as atom I. If all atoms have 3 substituents, this would mean 6 bonds + 6 centers = 12 improper torsions. Is that the minimum number of improper torsions required to keep a 6-membered ring flat? What about a 6-membered ring represented with united atoms? Thanks for reading so far (and many thanks for replying!). Regards, Pablo -- Pablo Englebienne, PhD Institute of Complex Molecular Systems (ICMS) Eindhoven University of Technology, TU/e PO Box 513, HG -1.26 5600 MB Eindhoven, The Netherlands Tel +31 40 247 5349 ___ gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users Please search the archive at http://www.gromacs.org/search before posting! Please don't post (un)subscribe requests to the list. Use the www interface or send it to gmx-users-requ...@gromacs.org. Can't post? Read http://www.gromacs.org/mailing_lists/users.php
