[gmx-users] further discussion on the mdrun -append function
Dear Acoot: You should be able to answer this one yourself. Moreover, you are doing yourself a disservice by relying on the mailing list to do your work for you because you will eventually need to learn how to find answers to these things on your own. Please remember the following: 1. use a title that matches your question, starting a new thread if you have a new question 2. it may sound harsh, but questions like this one in which you have obviously not even tried to find the answer yourself for more than a couple of minutes tend to annoy the very people that you may want help from later on with other issues. 3. whenever you post, please show how you have tried to solve the problem yourself. You may find that in the process of writing such a question you end up solving the problem yourself. 4. read the manual. Chris. -- original message -- Thanks for the detailed explaination. Will you please explain the function of -n index.ndx in grompp -f md.mdp -n index.ndx -p topol.top -c minimized.gro -o md-1ns.tpr with some specific examples? Cheers, Acoot -- 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] Different results from identical tpr after MD
Dear Acoot: The idea of convergence is this: start a large number of simulations from different conformations, analyze some quantity over time in each simulation, and when the deviation of the average value of that quantity from each separate simulation is less than the time-variance within individual simulations, then you can imply that the simulations have converged -- that is, the results of independent simulations which started as different are now similar. There is a huge body of work that uses a single simulations and evaluates its so-called convergence using some assumptions and special methods. That can also be very useful, but I find it informative to think of convergence in its standard non-scientific dictionary definition as the coming together of previously disparate things. For simulations, my working definition is this: a set of simulations has converged the value of some variable when the simulations were initiated from sufficiently distinct conformational basins and then, over time, the ensemble distribution of the time-averages of the specified variable has a variance that is the same as the mean time-averaged variation within independent simulations. The weak point here is the part about sufficiently distinct conformations, but I am not sure that this can be stated less vaguely in the general case. Chris. -- original message -- Hi Justin, Can you give me your definition of converged MD and unconverged MD? Cheers, Acoot -- 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] How to compile the template.c when GMX is builded with cmake
Dear Users: I was previously able to compile template.c after I compiled gromacs with autoconf. I was unable to compile templae.c, however, I used cmake to compile gromacs. This is from gromacs-4.5.4. I tried cmake . in the template directory with no success: cmake .gpc-f102n084-$ cmake . -- The C compiler identification is Intel -- The CXX compiler identification is Intel -- Check for working C compiler: /home/scinet/gpc/intel/ics/composer_xe_2011_sp1.6.233/bin/intel64/icc -- Check for working C compiler: /home/scinet/gpc/intel/ics/composer_xe_2011_sp1.6.233/bin/intel64/icc -- works -- Detecting C compiler ABI info -- Detecting C compiler ABI info - done -- Check for working CXX compiler: /home/scinet/gpc/intel/ics/composer_xe_2011_sp1.6.233/bin/intel64/icpc -- Check for working CXX compiler: /home/scinet/gpc/intel/ics/composer_xe_2011_sp1.6.233/bin/intel64/icpc -- works -- Detecting CXX compiler ABI info -- Detecting CXX compiler ABI info - done -- checking for module 'libmd' -- package 'libmd' not found CMake Error at CMakeLists.txt:42 (message): libmd not found, source GMXRC. -- Configuring incomplete, errors occurred! ### For this test, I used the template.c that came with gromacs. This issue has previously been posted to the list (see below) but I was unable to find any answers to that previous querry. I also looked at the README that came in the template directory and find it to be apparently autoconf specific, with no directions for cmake. Thank you, Chris. quoting http://www.mail-archive.com/gmx-users@gromacs.org/msg39347.html How to compile the template.c when GMX is builded with cmake Bert Thu, 31 Mar 2011 10:15:22 -0700 Dear all, When GMX is builded with cmake, how to compile the template.c? I used to make the template.c by the command make -f MakefileXXX, but it can not work in version 4.5. Thanks for your suggestions. Best regards, Bert -- -- 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] on the Umbrella Sampling on-line tutorial
Dear Acoot: I'll reply to general topics, not to the tutorial in particular. If the opening is not large enough to allow the peptide to exit, then the surrounding protein will need to change its conformation to permit unbinding. This can happen, but you need (a) to have sufficiently long sampling times to permit the required conformational change and (b) starting structures in which the peptide is near the umbrella center and do not crash during simulation. The US method is always valid, but it is not always the best choice. You might try using the free energy code to implement a thermodynamic cycle. Nevertheless, one imagines that the protein does actually need to open up during peptide binding and so you will still need to sample protein opening/closing in order to obtain equilibrium (i.e. correct) binding free energies because you need to sample the unbound protein at equilibrium. That is to say that a thermodynamic cycle may appear converged when in fact it is not, because you have not converged the unbound state of the protein. In any event, be careful with your convergence analysis. This would be a good system for which to attempt both US and double-decoupling approaches and use the results of each to ensure that you are not missing some important conformational states. That said, I highly doubt that it is possible to converge the free energies of induced-fit peptide-protein binding with any available atomistic computational method using contemporary computational resources. The lower bound of required sampling times is certainly on the order of 10 us per umbrella and I bet that the actual value is a few orders of magnitude larger. I have less experience with the free energy code than I do with US, but I suspect that the required sampling times are also very long for a system like this. Chris. -- original message -- Dear All, I planned to use the method introduced in the Umbrella Sampling on-line tutorial of Justin Lemkul (http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/index.html). But if a peptide is surrounded by a protein, which means the opening of the protein-complex is not large enough to allow the peptide to leave the protein without significantly breaking the conformation of the protein in the protein-peptide complex, is the Umbrella Sampling method still valid for the binding energy calculation? Will you please also show me in which part of the tutorial the direction of pull-apart is defined? We should process it in a direction the peptide can leave the protein, not the direction protein will bind the peptide much strongly. I am looking forward to getting a reply from you. Cheers, Acoot -- 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] Umbrella sampling along Radius of gyration
It is not clear to me how one would do this with MD. The only thing that I can think of doing in gromacs is to create a virtual particle that is placed at the center of the protein and then apply forces along the vector from this virtual atom to each of the Ca atoms. You would need to modify the gromacs source code so that the Rg was calculated at each step and the magnitude of each force is scaled appropriately such that you get a harmonic potential about the desired value of Rg. It should be easier to do with MC (although that's a ways off for gromacs unless I've missed something). Some people appear to have done exactly what you want with MD. I presume that they used charmm. http://www.pnas.org/content/94/19/10161.long (and their reference 22). Chris. -- original message -- Hi, Is there any way in gromacs to use radius of gyration of a polymer as reaction coordinate for umbrella sampling ? Thanks Sanku -- 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] Not able to continue with Equilibration
I didn't follow this whole thread, but I sometimes need to turn off all constraints when doing minimization. In fact, for that reason I entirely stopped ever using restraints during energy minimization. In extreem cases, I have had success also by forcing the water to be flexible with a -DFLEXIBLE (or whatever is appropriate for your water model; check the .itp) statement in the .mdp file. Once EM is done, use rigid water and restraints and everything works out ok. Chris. -- original message -- Hallo Felix, thank you for your answer. I tried the constraints = h-bonds but no change in the output. If I look at the step.pdb file that is produced after the running I have some strange outcome. For example some of my atoms are not recognized as part of my protein any more and my structure is destroied. Am I using some strange parameters for nvt without realizing it? I've started from the mdp file provided by the lysozyme tutorial non the Gromacs webpage. if anyone has any ideas it is more than welcome. Francesca integrator = md; leap-frog integrator nsteps = 1 ; 2 * 1 = 20 ps dt = 0.002 ; 2 fs ; Output control nstxout = 1 ; save coordinates every 0.2 ps nstvout = 1 ; save velocities every 0.2 ps nstenergy = 50; save energies every 0.2 ps nstlog = 50; update log file every 0.2 ps ; Bond parameters unconstrained_start = no ;continuation = no; first dynamics run constraint_algorithm = lincs; holonomic constraints constraints = h-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 = 3.0 ; short-range neighborlist cutoff (in nm) rcoulomb= 3.0 ; short-range electrostatic cutoff (in nm) rvdw= 3.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 = berendsen ; 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 = 300 300 ; reference temperature, one for each group, in K ; Pressure coupling is off pcoupl = no; no pressure coupling in NVT ; Periodic boundary conditions pbc = xyz ; 3-D PBC ; Dispersion correction DispCorr= EnerPres ; account for cut-off vdW scheme ; Velocity generation gen_vel = yes ; assign velocities from Maxwell distribution gen_temp= 300 ; temperature for Maxwell distribution gen_seed= -1; generate a random seed 2012/3/29 Rausch, Felix frau...@ipb-halle.de: [Hide Quoted Text] Hello again, Well, I once had problems with simulations crashing randomly during production runs (sometimes after tens of nanoseconds) with the LINCS warnings you described. Switching LINCS from all-bonds to only h-bonds did the trick for me, although I never exactly figured out why. Maybe it's worth a try for you, too? Cheers, Felix -Ursprüngliche Nachricht- Von: gmx-users-boun...@gromacs.org [mailto:gmx-users-boun...@gromacs.org] Im Auftrag von francesca vitalini Gesendet: Donnerstag, 29. März 2012 15:15 An: Discussion list for GROMACS users Betreff: Re: [gmx-users] Not able to continue with Equilibration Hi! I'm having a similar problem. I have a dimer solvated in a big box of water plus ions that I have managed to minimize correctly (see output of minimization at the end) but when I try to run NVT equilibration (see later) I get LINCS warnings(see below) refearred to atoms which are not in a cluster or in a strange position. I have added dihedral restraints on them but still the same type of error. I'm using GROMACS 3.3.1. I have tried to switch to a newer version of GROMACS but still the same error. Does anyone have any suggestions? Thanks Francesca MINIMIZATION OUTPUT Steepest Descents converged to Fmax 1000 in 681 steps Potential Energy = -1.9597512e+07 Maximum force = 2.4159846e+02 on atom 13087 Norm of force = 2.1520395e+04 MINIMIZATION.MDP define = -DEFLEXIBLE ; flexible water integrator = steep ; Algorithm (steep = steepest descent minimization) emtol = 1000.0 ; Stop minimization when the maximum force 1000.0 kJ/mol/nm emstep = 0.001 ; Energy step size nsteps = 5 ; Maximum number of (minimization) steps
[gmx-users] crazy temperatures
I disagree. What one is generally trying to obtain with elevated temperatures is enhanced sampling, not temperature-dependent properties. I believe that even TIP4P-EW is not very good at getting the properties of water correct at 600 K, temperatures that are commonly used during replica exchange simulations (not to mention that nobody has any idea how accurate protein forcefields are at temperatures other than the one at which they were parameterized). So I think that doing simulations at massively elevated temperatures can possibly be useful. That said, while doing simulated annealing, I have found previously using charmm that once you get to about 3,000 K you will get chiral inversions that can not resolve at lower temperature. This is because our improper dihedral terms only maintain the given chirality, rather than favouring one over the other. To address your question directly, I believe that chiral inversions will be a big problem for you at 20,000 K. Obviously you also have simulation stability issues, but one presumes that you could resolve those by using a small enough timestep. Chris. -- original message -- At that temperature most matter is going to be a plasma, not many bonds to be simulated and a lot of free electrons. Warren Gallin On 2012-03-28, at 4:43 PM, Mark Abraham wrote: [Hide Quoted Text] On 29/03/2012 9:39 AM, Asaf Farhi wrote: Dear GMCS users Hi. Does anyone know if MD at 2K is feasible? Please start new email threads rather than hijacking old ones. I doubt anybody knows the answer to your question. Force fields are parameterized to reproduce data at around 300K. I can't imagine any possible use for simulating an MM force field at a temperature hotter than the sun. Mark -- -- 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] Announcement: Large biomolecule benchmark report
You should absolutely publish this. it would be of great interest. You can mitigate your chances of running into problems with the overview by sending a version of the manuscript to the developers of each software and asking them to provide a short paragraph, each of which you could include in a final section of responses from the developers. A manuscript such as this (and indeed the information that you have already made available) will be very useful for many reasons. One reason is that when new PhD students start learning about simulations, they tend to use the package that has been adopted by their research group and trust the (probably biased and partially uninformed) statements of their senior colleagues. Chris. -- original message -- Thanks a a lot to you and also to Szilárd for your feedback and encouragement. I am very happy to see that this work is indeed useful especially to developers. We have no plans to make this into a 'proper' publication. I am not sure how much interested the simulation community would be because, to be honest, I have no overview what has been done in this area (besides the few benchmarks studies I have cited). Thanks again, Hannes. On Thu, 15 Mar 2012 22:02:21 +0100 David van der Spoel sp...@xray.bmc.uu.se wrote: [Hide Quoted Text] On 2012-03-15 14:37, Hannes Loeffler wrote: Dear all, we proudly announce our third benchmarking report on (large) biomolecular systems carried out on various HPC platforms. We have expanded our repertoire to five MD codes (AMBER, CHARMM, GROMACS, LAMMPS and NAMD) and to five protein and protein-membrane systems ranging from 20 thousand to 3 million atoms. Please find the report on http://www.stfc.ac.uk/CSE/randd/cbg/Benchmark/25241.aspx where we also offer the raw runtime data. We also plan to release the complete joint benchmark suite at a later date (as soon as we have access to a web server with sufficient storage space). We are open to any questions or comments related to our reports. It looks very interesting, and having benchmarks done by independent researchers is the best way to avoid bias. The differences are quite revealing, but it's also good that you point to problems compiling gromacs. Is this going to be submitted for publication somewhere too? Thanks for doing this, it must have been quite a job! Kind regards, Hannes Loeffler STFC Daresbury -- Scanned by iCritical. -- 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] intel grompp with pathscale mdrun
Dear users: can I use a .tpr file created with an intel icc compilation of grompp and then do mdrun under a pathscale compilation of mdrun ? (same version of gromacs) I'm wondering if there would be some strange behaviour. It seems like it should be ok, but I wanted to be sure. I ask because I can only compile mdrun with the pathscale compiler because make install hangs while make install-mdrun works fine. make install is hanging on making g_rms.o as previously posted: http://lists.gromacs.org/pipermail/gmx-users/2011-December/066682.html . I tried this again with the 4.0.12.1 version of the compiler and have the same result. Thank you, Chris. -- 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] intel grompp with pathscale mdrun
Thank you Matthew. your suggestion for the hanging .o compilation sounds good. I didn't try it because in the end the intel compiler produced the fastest executables in any event. For those interested, here are my benchmarking speeds for one of my simulation systems (270,000 atoms) on 48 cores of an AMD magny-cours cluster at 2.1 GHz. I have no idea why the intel icc compiler outperforms the pathscale and pgi compilers on AMD computers, but that seems to be the case. intel icc 12.0.5.220 = 3.733 ns/day intel icc 12.0.5.220 with -msse3 = 3.765 ns/day pathscale 4.0.12.1 = 3.467 ns/day pgi 11.8 = 3.092 ns/day pgi 11.8 with -tp istanbul -fast = 3.156 ns/day Again, thank you, Chris. -- original message -- As far as compilation hanging...maybe hand-compile that .o with less aggressive optimization flags, then try make again? MZ On Tue, Dec 6, 2011 at 2:20 PM, Chris Neale chris.neale at utoronto.ca wrote: Dear users: can I use a .tpr file created with an intel icc compilation of grompp and then do mdrun under a pathscale compilation of mdrun ? (same version of gromacs) I'm wondering if there would be some strange behaviour. It seems like it should be ok, but I wanted to be sure. I ask because I can only compile mdrun with the pathscale compiler because make install hangs while make install-mdrun works fine. make install is hanging on making g_rms.o as previously posted: http://lists.gromacs.org/pipermail/gmx-users/2011-December/066682.html . I tried this again with the 4.0.12.1 version of the compiler and have the same result. Thank you, Chris. -- gmx-users mailing listgmx-users at 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-request at gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists Previous message: [gmx-users] intel grompp with pathscale mdrun Next message: [gmx-users] dihedral format in top file? Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] -- 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] Half double pair list method in GROMACS [update]
What you have done seems alright. I didn't look closely enough to be sure though. One of the good things about this method is that you can easily test it yourself. To do this, create two different .gro files, one containing the atoms from one ff and the other containing the other atoms. For each, do separate zero-step mdrun evaluations of their energies under (a) their original ff, and (b) the combined HEDP ff that you have constructed. The energy evaluations should be the same once you account for rounding errors. Then all that is left is for you to be certain that you didn't cause any problems for nonbonded interactions. Note that you'll need to return to the original atomtypes for the first half of this test. Note that the only thing that the HEDP method is intended to perturb is the 1-4 interactions (and by perturb I mean that they will now be correct). Chris. -- original message -- Hi Stephane Chris, I followed all the threads posted by you two. I have a protein using ff99SB and GLYCAM for sugars. I have a disaccharide bound to protein. In xleap of AMBERTOOLS, I use the GLYCAM and ff99SB to generate the topology and coordinate files. I did the tests as Chris' original posts and by Stephane. The 1-4 interaction terms match for sugar alone and protein alone with HEDP method. When I generate topology and coordinate files with xleap of AMBER, there are three atom types that are common to protein and sugar. The atom types for my case are H1, OH and HO. Since for protein pairtypes using ff99SB, the epsilon has to be divided by 10 and pairs section be replicated five times and for sugar the epsilon in the pairtypes be divided by six and replicated six times, I am a little concerned about the three atomtypes that are common. So what I did was to change the atomtypes of sugar as H1S, OHS HOS for sugar and H1, OH and HO for protein and I made the changes accordingly in the pairtypes section for protein and sugar. Is this a valid approach? Any suggestion will be helpful. To make things little more clear: H1H10. 0. A 2.47135e-01 6.56888e-02 ;originally obtained using amb2gmx.pl H1S H1S 0. 0. A 2.47135e-01 6.56888e-02 ;Glycam Hydrogen of Sugar ( I changed this so that the common atom types be separated) In the ffnonbonded_complex_mod.itp: ;;using combination rule of 2 [ pairtypes ] ;;for protein H1 H1 1 0.247135000 0.006568880 ;the epsilon is divided by 10 ;;for sugar H1S H1S 1 0.24713500 0.010948133 ;the epsilon is divided by 6 Thanks for your time, Regards Sai On Mon, Sep 5, 2011 at 11:33 AM, ABEL Stephane 175950 Stephane.ABEL at cea.frwrote: Dear All, Below a little update and results about the application of half double pair list method to scale properly the Coulombic 1-4 interactions in case of a system where the AMBER99SB (fudgeLJ=0.5 and fudgeLJ=0.8333) and GLYCAM06 (fudgeLJ=1.0 and fudgeLJ=1.0) force fields are combined. I have followed the 4 steps described in [1] and used the following values in my forcefield.itp file [ defaults ] ; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ 1 2 yes 1.0 0.16 #include ffnonbonded_mod.itp ;#include ffnonbonded.itp #include ffbonded.itp I used two different topology files for the glycolipid (bDM) and the peptide. One with (*_mod.itp) with pair list parameters duplicated 6 times (bDM) and 5 times (peptide) and with single pair list (*_no_mod.itp) as decribed in [1]. TESTING: Three 3 different systems were examined: A. A first system containing 1 glycolipid (bDM) in water cubic box B. A second system with 1 peptide in TIP3P water C. And a third system with 1 peptide and 1 glycolipid in water cubic box To obtain the glycolipid and peptide energy pairs, I did one step of MD in NVT ensemble with the *.mdp file given in [2] with different energygrps and tc_grps. For 1. energygrps and tc_grps = bDM SOL For 2. energygrps and tc_grps = Protein SOL For 3. energygrps and = Protein bDM SOL bDM/water system Test_A1 ## Control with GLYCAM force field fudgeLJ fudgeQQ parameters and the *_no_mod.itp file : ##; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ ## 1 2 yes 1.0 1.0 Epot (kJ/mol)Coul-SR LJ-SRCoul-14 LJ-14 bDM-bDM -4.00855e+02 -3.71401e+012.03406e+032.79234e+02 Test_A2 with the topology *_mod.itp file and the directive ##; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ ## 1 2 yes 1.0 0.16 Epot (kJ/mol)Coul-SR LJ-SRCoul-14 LJ-14 bDM-bDM -4.00855e+02 -3.71401e+012.03406e+032.79234e+02 Test_A3 with the topology *_no_mod.itp file and the directive ##; nbfunccomb-rule gen-pairs
[gmx-users] Failed to lock .log. Already running simulation?
Dear Users: I have 50 simulations that are all the same, except with different random seeds for velocities. All were running fine for 24 hours. I canceled the running jobs and resubmitted them as part of beta testing a new cluster. All 50 started. I then canceled one of these jobs soon after starting it and then started it again pretty quickly (possibly too quickly). This restart now gave me the error: Fatal error: Failed to lock: continue.log. Already running simulation? For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors I found this post about this possibly being related to the Lustre filesystem: http://lists.gromacs.org/pipermail/gmx-users/2010-November/056173.html But I am not sure how to figure out if that is being used. Here is the output from mount: [nealechr@ip13-mp2 50]$ mount /dev/mapper/hddvg-root on / type ext4 (rw) proc on /proc type proc (rw) sysfs on /sys type sysfs (rw) devpts on /dev/pts type devpts (rw,gid=5,mode=620) tmpfs on /dev/shm type tmpfs (rw) /dev/md0 on /boot type ext4 (rw) /dev/mapper/hddvg-home on /home type ext4 (rw,usrquota,grpquota) /dev/md2 on /ltmp type ext4 (rw) /dev/mapper/hddvg-opt on /opt type ext4 (rw) none on /ramdisk type tmpfs (rw,nosuid,nodev) none on /var/tmp type tmpfs (rw,noexec,nosuid,nodev,size=10) none on /proc/sys/fs/binfmt_misc type binfmt_misc (rw) none on /ipathfs type ipathfs (rw) sunrpc on /var/lib/nfs/rpc_pipefs type rpc_pipefs (rw) nfsd on /proc/fs/nfsd type nfsd (rw) none on /tmp type tmpfs (rw,noexec,nosuid,nodev,size=10) 10.4.215.201@o2ib:/lustre01 on /mnt/scratch01 type lustre (rw,_netdev,flock) Also, it seems unlikely to be system related because the other 49 runs are going just fine. I did a ls -la to see if there was some hidden file to indicate the lock but could not find any (I have no idea how such a lock would work or be detected). I deleted the .log file, but then I get the error: Fatal error: File appending requested, but only 3 of the 4 output files are present Moving everything to a new directory and then copying it back (including the original .log file) allowed me to run the simulation. Did I do something incorrectly, or is this a bona-fide problem? Thank you, Chris. -- 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] compiling with the PGI compiler
I am using a new cluster of Xeons and, to get the most efficient
compilation, I have compiled gromacs-4.5.4 separately with the intel,
pathscale, and pgi compilers.
With the pgi compiler, I am most concerned about this floating point
overflow warning:
...
[ 19%] Building C object src/gmxlib/CMakeFiles/gmx.dir/trajana/trajana.c.o
[ 19%] Building C object
src/gmxlib/CMakeFiles/gmx.dir/trajana/centerofmass.c.o
[ 19%] Building C object src/gmxlib/CMakeFiles/gmx.dir/trajana/nbsearch.c.o
PGC-W-0129-Floating point overflow. Check constants and constant
expressions
(/home/nealechr/exe/pgi/gromacs-4.5.4/source/src/gmxlib/trajana/nbsearch.c:
166)
PGC/x86-64 Linux 11.8-0: compilation completed with warnings
[ 19%] Building C object
src/gmxlib/CMakeFiles/gmx.dir/trajana/displacement.c.o
[ 19%] Building C object src/gmxlib/CMakeFiles/gmx.dir/trajana/position.c.o
...
There are also a lot of warnings about a type cast, which seems less
like a real problem:
...
[ 85%] Building C object src/mdlib/CMakeFiles/md.dir/gmx_qhop_xml.c.o
[ 85%] Building C object src/mdlib/CMakeFiles/md.dir/nsgrid.c.o
[ 85%] Building C object src/mdlib/CMakeFiles/md.dir/stat.c.o
[ 86%] Building C object src/mdlib/CMakeFiles/md.dir/mdebin_bar.c.o
[ 86%] Building C object src/mdlib/CMakeFiles/md.dir/genborn_sse2_double.c.o
[ 86%] Building C object src/mdlib/CMakeFiles/md.dir/pull.c.o
PGC-W-0095-Type cast required for this conversion
(/home/nealechr/exe/pgi/gromacs-4.5.4/source/src/mdlib/pull.c: 1213)
PGC-W-0095-Type cast required for this conversion
(/home/nealechr/exe/pgi/gromacs-4.5.4/source/src/mdlib/pull.c: 1213)
PGC/x86-64 Linux 11.8-0: compilation completed with warnings
[ 86%] Building C object src/mdlib/CMakeFiles/md.dir/force.c.o
PGC-W-0095-Type cast required for this conversion
(/home/nealechr/exe/pgi/gromacs-4.5.4/source/src/mdlib/force.c: 519)
PGC/x86-64 Linux 11.8-0: compilation completed with warnings
[ 86%] Building C object src/mdlib/CMakeFiles/md.dir/fft5d.c.o
[ 86%] Building C object src/mdlib/CMakeFiles/md.dir/tables.c.o
[ 86%] Building C object src/mdlib/CMakeFiles/md.dir/gmx_parallel_3dfft.c.o
[ 86%] Building C object src/mdlib/CMakeFiles/md.dir/qmmm.c.o
[ 88%] Building C object src/mdlib/CMakeFiles/md.dir/vcm.c.o
PGC-W-0095-Type cast required for this conversion
(/home/nealechr/exe/pgi/gromacs-4.5.4/source/src/mdlib/vcm.c: 83)
PGC/x86-64 Linux 11.8-0: compilation completed with warnings
[ 88%] Building C object src/mdlib/CMakeFiles/md.dir/sim_util.c.o
...
###
I compiled like this:
module purge
module load pgi64/11.8 openmpi_pgi64/1.4.3_ofed
module load cmake/2.8.5
export CCDIR=/opt/pgi/linux86-64/11.8/bin
## set the location of the single precision FFTW isntallation
export FFTW_LOCATION=/home/nealechr/exe/pgi/fftw-3.1.2/exec
# Nothing below this line usually needs to be changed
export CXX=pgCC
export CC=pgcc
cmake ../source/ \
-DFFTW3F_INCLUDE_DIR=$FFTW_LOCATION/include \
-DFFTW3F_LIBRARIES=$FFTW_LOCATION/lib/libfftw3f.a \
-DCMAKE_INSTALL_PREFIX=$(pwd) \
-DGMX_X11=OFF \
-DCMAKE_CXX_COMPILER=${CCDIR}/pgCC \
-DCMAKE_C_COMPILER=${CCDIR}/pgcc \
-DGMX_PREFER_STATIC_LIBS=ON \
-DGMX_MPI=OFF
make
make install
##
I don't get any similar errors with the intel or pathscale compilers
(although intel gives me icc: command line warning #10159: invalid
argument for option '-m' a lot) and the pathscale compilation appears
to be hung on [ 84%] Building C object
src/tools/CMakeFiles/gmxana.dir/gmx_rms.c.o, but I'll post a separate
ticket for that if it remains hung for a long time.
If anybody has run into this warning, or knows enough to be sure that I
don't need to worry about it during mdrun (it seems to be in an analysis
file, but I'm not entirely sure), then I would be happy to hear about
it. A gmx-users search for PGI returned zero results. I saw something
here, but it was not very specific about the problem with pgi (
http://www.levlafayette.com/node/175 ).
Thank you,
Chris.
--
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] diffusion of the water at the micelle surface
Well, the positive way of looking at this is that it appears that nobody has ever done it before. If somebody has done it (in charmm for instance) then you might be able to convert the format of your .xtc and use their tools to analyze it. Good luck, Chris. -- original message -- Thank you very much for your reply. Indeed I have read several papers where the diffusion of water at the membrane surface have been computed. Since the diffusion of the interfacial water is an useful properties to examine the micelle/surface irregularities, I would hope that this option exist in gromacs. Unfortunately, it is not the case, so i will try your suggestions . Stephane -- Message: 6 Date: Sat, 26 Nov 2011 10:17:10 -0500 From: chris.neale at utoronto.ca Subject: [gmx-users] diffusion of the water at the micelle surface To: gmx-users at gromacs.org Message-ID: 2026101710.bbl5wl376s04ccww at webmail.utoronto.ca Content-Type: text/plain; charset=ISO-8859-1; DelSp=Yes; format=flowed When people do this for lipid bilayers, they compute depth-dependent diffusion profiles (often diffusion is computed separately for lateral diffusion and diffusion along the bilayer normal). Sounds like you might do something similar. I doubt that the standard gromacs tools will do this for you. If you don't hear from anybody about how to do this, then I'd suggest that you simply use g_dist to get the time-dependent distance for each water molecule and then use g_traj to output the coordinates of each water molecule and then script it yourself after reading one of the papers where people compute depth-dependent diffusion profiles for a lipid bilayer. Chris. -- original message -- I would like to compute the translational diffusion around the micelle surface. I know that I can select the water molecules at x distance of the micelle surface with g_select (right ?) but how to use this file generated by g_select to compute de diffusion, since the index and/or the number of water will change with the simulation time . Thank you for your response -- 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] calculating the pmf in constrained force simulations
Do you mean that you do constrained position simulations and want to know how to process the force? If so, read about thermodynamic integration (TI). I mostly work with US and position restraints, but for absolute constraints I believe that you should take the average force at each constrained position and do a trapezoidal integration of the mean force. Surely somebody has written a paper on this that you can read. Chris. -- original message -- Hi I do constrained force simulations and i have the pullf.xvg and pullx.xvg files.. I want to know how to calculate the force for each simulation and i how to integrate the forces.. can some one help me. -- 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] to gro or not to gro
It's more useful if you provide more information. What was the .pdb file (can I download it from the pdb databank?) was there water? what version of gromacs? was it compiled in double or single precision? what were your mdp parameters? -- original message -- There is something not quite right here. I grompp a gro file created from a pdb file using pdb2gmx and run MD to do a single point energy calculation. I then grompp the original pdb file and run MD to do a single point energy calculation. The values of potential energy values are different by about 3 kJ/mol for a small molecule. These things should carry a bad for your health sign :D -- 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] diffusion of the water at the micelle surface
When people do this for lipid bilayers, they compute depth-dependent diffusion profiles (often diffusion is computed separately for lateral diffusion and diffusion along the bilayer normal). Sounds like you might do something similar. I doubt that the standard gromacs tools will do this for you. If you don't hear from anybody about how to do this, then I'd suggest that you simply use g_dist to get the time-dependent distance for each water molecule and then use g_traj to output the coordinates of each water molecule and then script it yourself after reading one of the papers where people compute depth-dependent diffusion profiles for a lipid bilayer. Chris. -- original message -- I would like to compute the translational diffusion around the micelle surface. I know that I can select the water molecules at x distance of the micelle surface with g_select (right ?) but how to use this file generated by g_select to compute de diffusion, since the index and/or the number of water will change with the simulation time . Thank you for your response Stephane -- 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] to gro or not to gro
There is rounding because of angstroms vs nanometers and both files maintain 3 decimal places (see below). The intention of pdb2gmx is not to get a .gro , but to get a .top or .itp file. I see your point, but I don't think it matters. If you really need that precision, then I'd think that you'd be using double precision gromacs (which would still round the output .gro from pdb2gmx and not solve this issue... just saying). gpc-f101n084-$ tail a.gro 2GLYHA2 11 -0.426 0.215 0.193 2GLY C 12 -0.446 0.390 0.072 2GLY O 13 -0.368 0.440 -0.013 3NME N 14 -0.532 0.478 0.154 3NME H 15 -0.618 0.488 0.104 3NMECH3 16 -0.468 0.607 0.173 3NME HH31 17 -0.375 0.594 0.227 3NME HH32 18 -0.533 0.670 0.234 3NME HH33 19 -0.447 0.661 0.081 0.47330 0.80470 0.29640 gpc-f101n084-$ tail a.pdb ATOM 11 HA2 GLY 2 -4.263 2.147 1.932 ATOM 12 C GLY 2 -4.461 3.901 0.720 ATOM 13 O GLY 2 -3.677 4.400 -0.128 HETATM 14 N NME 3 -5.316 4.776 1.535 HETATM 15 H NME 3 -6.182 4.876 1.044 HETATM 16 CH3 NME 3 -4.685 6.072 1.731 HETATM 17 1HH3 NME 3 -3.754 5.941 2.272 HETATM 18 2HH3 NME 3 -5.328 6.697 2.341 HETATM 19 3HH3 NME 3 -4.465 6.615 0.809 END -- original message -- There was no water. Version 4.5.5. Single precision. This is what I ran: pdb2gmx -f 1.pdb -p g1 -o g1 grompp -f 1.mdp -c g1.gro -o g1.tpr -p g1.top mdrun -v -s g1.tpr -c ag1.gro -g g1.log -e g1.ene Potential energy: 1.19780e+02 kJ/mol grompp -f 1.mdp -c 1.pdb -o g1.tpr -p g1.top mdrun -v -s g1.tpr -c ag1.gro -g g1.log -e g1.ene Potential energy: 1.15997e+02 kJ/mol The mdp file: integrator = md nsteps = 0 nstcomm = 1 nstxout = 1 nstvout = 1 nstfout = 0 nstlog = 1 nstenergy= 1 nstxtcout= 0 nstlist = 0 ns_type = simple pbc = no rlist= 0 rcoulomb = 0 rvdw = 0 gen_vel = no unconstrained-start = yes lincs-warnangle = 30 The pdb file : HETATM1 1HH3 ACE 1 -1.449 0.109 -0.581 H HETATM2 CH3 ACE 1 -2.101 -0.268 0.202 C HETATM3 2HH3 ACE 1 -1.708 0.090 1.151 H HETATM4 3HH3 ACE 1 -2.109 -1.350 0.189 H HETATM5 C ACE 1 -3.491 0.269 0.000 C HETATM6 O ACE 1 -4.453 -0.401 -0.152 O ATOM 7 N GLY 2 -3.601 1.734 0.000 N ATOM 8 H GLY 2 -3.039 2.278 -0.623 H ATOM 9 CA GLY 2 -4.525 2.391 0.903 C ATOM 10 HA1 GLY 2 -5.539 2.051 0.696 H ATOM 11 HA2 GLY 2 -4.263 2.147 1.932 H ATOM 12 C GLY 2 -4.461 3.901 0.720 C ATOM 13 O GLY 2 -3.677 4.400 -0.128 O HETATM 14 N NME 3 -5.316 4.776 1.535 N HETATM 15 H NME 3 -6.182 4.876 1.044 H HETATM 16 CH3 NME 3 -4.685 6.072 1.731 C HETATM 17 1HH3 NME 3 -3.754 5.941 2.272 H HETATM 18 2HH3 NME 3 -5.328 6.697 2.341 H HETATM 19 3HH3 NME 3 -4.465 6.615 0.809 H END *In reply to chris.neale at utoronto.ca:* It's more useful if you provide more information. What was the .pdb file (can I download it from the pdb databank?) was there water? what version of gromacs? was it compiled in double or single precision? what were your mdp parameters? -- 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] to gro or not to gro
1. why repeat the calculations? If you're talking about simulations then there is no need to repeat them due to this. You will get different answers with the same starting coordinates if you simply change the initial velocities. If you're talking about instantaneous energy calculations then I suppose you might need to redo it, but they should be very quick, right? 2. The .gro files do not carry useless zeroes. you have it backwards... the gro files end up with fewer digits. 3. it's a little annoying to find out that you already knew the answer. Why not state that at the outset? Unless I misunderstand this point, this will mark the end of my comments since holding back information on purpose just wastes people's time. Chris. -- original message -- I already knew the reason. But I had to find this out hard way. Now facing a dreading prospect of repeating tons of calculations! Hence the request for the bad for the health sign. Btw, gromacs issues a lot of warnings, but not this one :D What a bright idea to switch from angstroms to nanometers! Now the gro files carry a lot of useless zeros. -- 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] Charged wall
You could create an atomistic charged wall by placing atoms in a plane or grid and using position restraints or freeze groups to keep them in place. You could them use the pull code to restrain part of the protein relative to an atom(s) in this atomistic wall. You could obtain your desired charge distribution by using a combination of helium and sodium or chloride atoms. Alternatively, you could create a new atom type with a fractional charge and whatever LJ parameters you want. Of course, you could also construct a more atomistically realistic wall.Some people have done similar things to study proteins absorbed onto electrodes: http://pubs.acs.org/doi/abs/10.1021/ja910707r http://www.sciencedirect.com/science/article/pii/S0013468609003119 Chris. -- original message -- Hi Gmx Users, I am wondering whether is it possible in Gromacs to build sth like a charged wall to which the last residue of my terminal protein will be attached and run the simulation? I mean that the protein cannot move through some border which is attached to and which is positively or negatively charged. Thank you Steven -- 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] umbrella sampling with pull=constraint
Dear Vijay: Can you please provide evidence for your claim that the harmonic potential is not applied properly, since you may decide to use pull=umbrella once you have set that up correctly. Importantly, movement out of the unit-cell is not a problem, as discussed a lot on this list. Nevertheless, you do need to worry about which images are used to derive the pulling forces. You can often do that by a judicious selection of pull_pbcatom (read about that on-list and in the manual). In some cases, however, pull_pbcatom can not assist enough and you are forced to make the box larger. None of this is alleviated by pull=constraint, which is why I'm not going to answer your question about your .mdp parameters at this point. Let's get the problems identified and solved first with pull=umbrella. Chris. -- original message -- Hello, I have done umbrella sampling with pull=umbrella and I found that the pulling group has high fluctuations and sometimes moving out of the periodic box. I think that the harmonic potential is not properly applied and thus the pulling group is not retained with the specified COM distance between the reference and pulling group. Can we use pull=constraint option to retain the pulling group within the COM distance? my pulling code is as bellow with restrain at 0.5nm distance from ref group. I just want to get some idea about this pull code modification. pull= constraint pull_geometry= distance pull_dim= N N Y pull_start = no pull_ngroups= 1 pull_group0 = Chain_A pull_group1 = Chain_B pull_init1 = 0.50 pull_rate1 = 0.0 pull_k1 = 1000 pull_nstxout= 500 pull_nstfout= 500 Regards, Vijay. -- 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] lipid bilayer and umbrella sampling
If you add position restraints to your DPPC molecules, then you are changing your effective order parameter. The PMF for the solute along the normal to a restrained bilayer will have a different shape than a PMF for the solute along the normal to an urestrained bilayer. Whether or not this is a problem will depend on the question that you are trying to answer. If the bilayer is falling apart, then you certainly need to do something. But if it is just locally ordering/disordering, then I don't see the big problem, besides the effect that this has on convergence (DOI: 10.1021/ct200316w). Chris. -- original message -- Dear Gromacs Users, I am trying to extract the potential of mean force of a small molecule in a DPPC bilayer. To this end, I applied the methodology described in an online manual written by Justin Lemkul. My problem is when I run biasing simulations of the molecule near the interface (DPPC/water), some lipid molecules move to the water phase. This has as a consequence a local disorder of the bilayer. Below is the parameters I employ for the pull code: pull = umbrella pull_geometry= distance pull_dim = N N Y pull_start = yes pull_ngroups = 1 pull_group0 = DPPC pull_group1 = DTC pull_rate1 = 0.0 pull_k1 = 1000 Is the position restraints on DPPC molecules a solution to my problem? Thanks in advance Best regards Giovani -- 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] source of opls Mg2+ parameters?
Thank you Tom! Aqvist_A12 = sqrt(gromacs_C12*10^12/4.184) Aqvist_A6 = sqrt(gromacs_C6*10^6/4.184) (equation verified based on the SPC water oxygen parameters, also listed in Aqvist). This reference (Aqvist, J (1990) J. Phys. Chem., 94 (21), pp 8021?8024) should probably be noted somewhere in ions.itp or in ffoplsaa.itp. That might help stop these non-opls parameters from being cited as opls parameters. Chris. -- original message -- Hi Chris, If I remember correctly the magnesium ion parameters are the Aquvist parameters (http://pubs.acs.org/doi/abs/10.1021/j100384a009). Cheers Tom chris.ne...@utoronto.ca wrote: [Hide Quoted Text] Dear users: does anybody know where the OPLS magnesium parameters are from? As far as I can tell, they are not in Jorgensen 1996 or Kaminski 2001, In spite of the fact that many simulation studies reference these papers for their magnesium opls parameters. In fact, I do not think that the opls mg2+ parameters in ions.itp are taken from any of the references listed in ffoplsaa.itp. I checked the following references with searches for magnesium, mg2+, and a general visual scan: ; W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, ; J. Am. Chem. Soc. 118, 11225-11236 (1996). ; W. L. Jorgensen and N. A. McDonald, Theochem 424, 145-155 (1998). ; W. L. Jorgensen and N. A. McDonald, J. Phys. Chem. B 102, 8049-8059 (1998). ; R. C. Rizzo and W. L. Jorgensen, J. Am. Chem. Soc. 121, 4827-4836 (1999). ; M. L. Price, D. Ostrovsky, and W. L. Jorgensen, J. Comp. Chem. (2001). ; E. K. Watkins and W. L. Jorgensen, J. Phys. Chem. A 105, 4118-4125 (2001). ; G. A. Kaminski, R.A. Friesner, J.Tirado-Rives and W.L. Jorgensen, J. Phys. Chem. B 105, 6474 (2001). From a web search, the TINKER program indicates that it uses OPLS Mg2+ parameters from an unpublished source: OPLS All-Atom Parameters for Organic Molecules, Ions, Peptides Nucleic Acids, July 2008 as provided by W. L. Jorgensen, Yale University during June 2009. These parameters are taken from those distributed with BOSS Version 4.8. (above quote from http://dasher.wustl.edu/ffe/distribution/params/oplsaa.prm ) although those tinker parameters (sigma=0.1450 and epsilon=3.9600) do not match the values for opls in the gromacs ions.itp (sigma=0.164447 and epsilon=3.66118). I did not check out BOSS directly (it costs money), but I did look at the manual on the Jorgensen web page, although that manual does not contain Mg2+ parameters. I also checked the gromos magnesium parameters that are distributed with gromacs and verified that these c6/c12 values do not match the sigma/epsilon opls values after a conversion with g_sigeps. Thank you very much for your time and assistance, Chris. -- Dr Thomas Piggot University of Southampton, UK. -- 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] source of opls Mg2+ parameters?
Dear users: does anybody know where the OPLS magnesium parameters are from? As far as I can tell, they are not in Jorgensen 1996 or Kaminski 2001, In spite of the fact that many simulation studies reference these papers for their magnesium opls parameters. In fact, I do not think that the opls mg2+ parameters in ions.itp are taken from any of the references listed in ffoplsaa.itp. I checked the following references with searches for magnesium, mg2+, and a general visual scan: ; W. L. Jorgensen, D. S. Maxwell, and J. Tirado-Rives, ; J. Am. Chem. Soc. 118, 11225-11236 (1996). ; W. L. Jorgensen and N. A. McDonald, Theochem 424, 145-155 (1998). ; W. L. Jorgensen and N. A. McDonald, J. Phys. Chem. B 102, 8049-8059 (1998). ; R. C. Rizzo and W. L. Jorgensen, J. Am. Chem. Soc. 121, 4827-4836 (1999). ; M. L. Price, D. Ostrovsky, and W. L. Jorgensen, J. Comp. Chem. (2001). ; E. K. Watkins and W. L. Jorgensen, J. Phys. Chem. A 105, 4118-4125 (2001). ; G. A. Kaminski, R.A. Friesner, J.Tirado-Rives and W.L. Jorgensen, J. Phys. Chem. B 105, 6474 (2001). From a web search, the TINKER program indicates that it uses OPLS Mg2+ parameters from an unpublished source: OPLS All-Atom Parameters for Organic Molecules, Ions, Peptides Nucleic Acids, July 2008 as provided by W. L. Jorgensen, Yale University during June 2009. These parameters are taken from those distributed with BOSS Version 4.8. (above quote from http://dasher.wustl.edu/ffe/distribution/params/oplsaa.prm ) although those tinker parameters (sigma=0.1450 and epsilon=3.9600) do not match the values for opls in the gromacs ions.itp (sigma=0.164447 and epsilon=3.66118). I did not check out BOSS directly (it costs money), but I did look at the manual on the Jorgensen web page, although that manual does not contain Mg2+ parameters. I also checked the gromos magnesium parameters that are distributed with gromacs and verified that these c6/c12 values do not match the sigma/epsilon opls values after a conversion with g_sigeps. Thank you very much for your time and assistance, Chris. -- 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] umbrella sampling convergence
The criteria are the same for any type of simulation. Generally, you must show that, as far as you can tell, the values that you derive are not going to change if you run the simulation a lot longer. Different quantities converge at different rates, so ideally you should check them all independently. In practice, people tend to pick one that they think will converge the most slowly and analyze the convergence of that. This could be your PMF. If it's not converged then you can either run longer or figure out what is relaxing slowly and include that degree of freedom in your reaction coordinate. Chris. -- original message -- What is the criteria for umbrella sampling convergence. If I am right, the pulling force and COM distance should be converged. I do not see force or COM distance convergence after 10ns of simulation. Regards, Vijay -- 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] Mismatching in Bergers DMPC bilayer lipid topology
Dear James: Next time, please specify exactly what you did in enough detail for somebody else to reproduce it, much as in a manuscript. e.g. there is no dmpc.gro in that website. I can guess what you did, but that is not ideal. I took a look at the files that I suppose you used and the atom names do appear to be mis-matched between the structure and the topology. That's is not necessarily a problem though. What happens if you set maxwarn to a large value during grompp and do a short mdrun? If that succeeds and looks normal, then I'd suggest to draw the structure of the lipid (about 50 atoms) and check out the topology by copying the names/charges/etc. onto your map. I, for one, have used the dmpc.itp from that website and found that it is ok. I accessed the file years ago, so there is no guarantee it's the same file, but my topology also has names like NTM. I never did use that DMPC structure file though. Chris. -- original message -- Is here anobody who worked with the dmpc bilayers obtained from http://moose.bio.ucalgary.ca/index.php?page=Structures_and_Topologies I've used this bilayers with the dmpc.itp and gromos96 ff with Berger lipids and during loading my bilayer.gro in the grompt I've obtained warning about mismatching in some atoms in .gro relatively .itp Warning: atom name 1 in topol_dmpc.top and dmpc.gro does not match (CN1 - CA) Warning: atom name 2 in topol_dmpc.top and dmpc.gro does not match (CN2 - CB) Warning: atom name 3 in topol_dmpc.top and dmpc.gro does not match (CN3 - CC) Warning: atom name 4 in topol_dmpc.top and dmpc.gro does not match (NTM - N) Warning: atom name 5 in topol_dmpc.top and dmpc.gro does not match (CA - CD) Warning: atom name 6 in topol_dmpc.top and dmpc.gro does not match (CB - CE) Warning: atom name 12 in topol_dmpc.top and dmpc.gro does not match (CC - CF) Warning: atom name 13 in topol_dmpc.top and dmpc.gro does not match (CD - CG) Warning: atom name 30 in topol_dmpc.top and dmpc.gro does not match (CE - CH) Warning: atom name 47 in topol_dmpc.top and dmpc.gro does not match (CN1 - CA) Warning: atom name 48 in topol_dmpc.top and dmpc.gro does not match (CN2 - CB) Warning: atom name 49 in topol_dmpc.top and dmpc.gro does not match (CN3 - CC) Warning: atom name 50 in topol_dmpc.top and dmpc.gro does not match (NTM - N) Warning: atom name 51 in topol_dmpc.top and dmpc.gro does not match (CA - CD) Warning: atom name 52 in topol_dmpc.top and dmpc.gro does not match (CB - CE) Warning: atom name 58 in topol_dmpc.top and dmpc.gro does not match (CC - CF) Warning: atom name 59 in topol_dmpc.top and dmpc.gro does not match (CD - CG) Warning: atom name 76 in topol_dmpc.top and dmpc.gro does not match (CE - CH) Warning: atom name 93 in topol_dmpc.top and dmpc.gro does not match (CN1 - CA) Warning: atom name 94 in topol_dmpc.top and dmpc.gro does not match (CN2 - CB) (more than 20 non-matching atom names) Does this warning harmfull ? :) JAmes - -- 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] umbrella sampling
You need to evaluate convergence yourself for any simulation. I suggest doing the whole thing twice (or more) with different starting conformations. Also, look at the PMF from block averaging (generate one PMF from the 0-2 ns data, another from the 2-4 ns data, and so on) and see if there is a systematic drift with time. In my opinion, if you really want to converge to within 5 kcal/mol without using additional restraints as Justin mentioned, then you should be expecting to do 100 ns per umbrella, and probably 1000 ns per umbrella. -- Thus additional restraints are a very good idea. To see about using additional restraints, check out this paper: THE JOURNAL OF CHEMICAL PHYSICS 125, 084902 (2006) and this one: Journal of Chemical Theory and Computation 3(4):1231-1235 (2007). They are not US, but the idea is the same. Chris. -- original message -- Hi Justin, Here is my complete code for umbrella sampling title = Umbrella pulling simulation define = -DPOSRES_2 integrator = md dt = 0.002 tinit = 0 nsteps = 250 ; 5 ns nstcomm = 10 nstxout = 5 ; every 100 ps nstvout = 5 nstfout = 5000 nstxtcout = 5000 ; every 10 ps nstenergy = 5000 constraint_algorithm= lincs constraints = all-bonds continuation= yes nstlist = 5 ns_type = grid rlist = 1.4 rcoulomb= 1.4 rvdw= 1.4 ; PME electrostatics parameters coulombtype = PME fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-5 optimize_fft= yes ; Berendsen temperature coupling is on in two groups Tcoupl = V-rescale tc_grps = Protein SOL tau_t = 0.5 0.5 ref_t = 300 300 ; Pressure coupling is on Pcoupl = Parrinello-Rahman pcoupltype = isotropic tau_p = 1.0 compressibility = 4.5e-5 ref_p = 1.0 ; Generate velocities is off gen_vel = no ; Periodic boundary conditions are on in all directions pbc = xyz ; Long-range dispersion correction DispCorr= EnerPres ; Pull code pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = Chain_B pull_group1 = Chain_A pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1000 pull_nstxout= 1000 pull_nstfout= 1000 I will do this umbrella sampling for 10ns, Do you think the COM distance can converge after long simulation time? ie back to the restraint position. Regards, Vijay. [Hide Quoted Text] Vijayaraj wrote: Hello, I am doing umbrella sampling to calculate PMF curve for the detachment of a terminal cyclic peptide with 8 aa's (CP) from the self-assembled cyclic peptide nanotube. I extracted the reaction coordinates staring from 5.5 ang COM distance between the terminal CP and the 2nd CP to 17.5 ang, I did umbrella sampling on 25 configurations (for 5ns) and the window size is 0.5 ang. I used the following pulling code for umbrella sampling, pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = Chain_B pull_group1 = Chain_A pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1000 pull_nstxout= 1000 pull_nstfout= 1000 I restrained the 2nd CP unit and the pull_rate1 is 0, so the COM distance between the chain_A (terminal) and chain_B (2nd CP) should be restrained. after 5ns of umbrella sampling, I calculated the COM distance between chain A and B, but it was not restrained, for the 5.5 ang COM distance restrain, the COM distance varies from 4.5 to 5.5 ang. and also the pulling cyclic peptide undergoes large conformational sampling. from the WHAM analysis I understood the sampling window is poorly represented. In addition to COM distance restrain, can I restrain the pulling CP's backbone atoms? so that the pulling groups large conformational sampling will be reduced. You could implement dihedral restraints to fix the backbone secondary structure, but I can't comment on the stability of trying to use these restraints in addition to the pull code. Seems like a lot going on at once, to me. Also consider the fact that 5 ns is an extremely short timeframe to gather meaningful data. Perhaps you just need more time in each window to equilibrate. At the shortest COM distance, your two molecules are still likely experiencing some interactions, and it may require a great deal of sampling in this window to converge the simulations. You haven't shown the rest of your .mdp file (always a good idea!), so we can only guess at whether or not your other settings should lead to a sensible result. -Justin -- 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
[gmx-users] RE:How to change the screen's output frequency?
Speaking of which, I think that all programs should by default list -hidden when they are run with -h so that users are aware that there are hidden options and know how to access them. The -h output could clearly indicate that such options are not fully tested, etc. Chris. Justin A. Lemkul jalemkul at vt.edu Mon Oct 3 15:56:45 CEST 2011 * Previous message: [gmx-users] RE:How to change the screen's output frequency? * Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] ?? wrote: On 2011-10-02 17:14, ?? wrote: How to change the screen's output frequency? ol 0.85 imb F 6% step 4188000, remaining runtime: 174 s vol 0.87 imb F 3% step 4188100, remaining runtime: 174 s vol 0.86 imb F 1% step 4188200, remaining runtime: 174 s vol 0.86 imb F 4% step 4188300, remaining runtime: 174 s vol 0.88 imb F 6% step 4188400, remaining runtime: 174 s vol 0.87 imb F 3% step 4188500, remaining runtime: 174 s vol 0.88 imb F 8% step 4188600, remaining runtime: 174 s vol 0.88 imb F 6% step 4188700, remaining runtime: 174 s vol 0.86 imb F 3% step 4188800, remaining runtime: 174 s vol 0.87 imb F 4% step 4188900, remaining runtime: 174 s vol 0.86 imb F 1% step 4189000, remaining runtime: 174 s vol 0.86 imb F 5% step 4189100, remaining runtime: 174 s vol 0.87 imb F 4% step 4189200, remaining runtime: 173 s vol 0.88 imb F 6% step 4189300, remaining runtime: 173 s vol 0.85 imb F 1% step 4189400, remaining mdrun -h I have done that ,but i did not fing the option! The option is hidden, try mdrun -h -hidden. You can use the -stepout option to set the output frequency. Or, you can turn off -v altogether. -Justin The following are the mdrun's options . Option Filename Type Description -s topol.tpr InputRun input file: tpr tpb tpa -o traj.trr Output Full precision trajectory: trr trj cpt -x traj.xtc Output, Opt. Compressed trajectory (portable xdr format) -cpi state.cpt Input, Opt. Checkpoint file -cpo state.cpt Output, Opt. Checkpoint file -cconfout.gro Output Structure file: gro g96 pdb etc. -e ener.edr Output Energy file -g md.log Output Log file -dhdl dhdl.xvg Output, Opt. xvgr/xmgr file -fieldfield.xvg Output, Opt. xvgr/xmgr file -tabletable.xvg Input, Opt. xvgr/xmgr file -tablep tablep.xvg Input, Opt. xvgr/xmgr file -tableb table.xvg Input, Opt. xvgr/xmgr file -rerunrerun.xtc Input, Opt. Trajectory: xtc trr trj gro g96 pdb cpt -tpitpi.xvg Output, Opt. xvgr/xmgr file -tpid tpidist.xvg Output, Opt. xvgr/xmgr file -eisam.edi Input, Opt. ED sampling input -eosam.edo Output, Opt. ED sampling output -j wham.gct Input, Opt. General coupling stuff -jobam.gct Output, Opt. General coupling stuff -ffout gct.xvg Output, Opt. xvgr/xmgr file -devout deviatie.xvg Output, Opt. xvgr/xmgr file -runav runaver.xvg Output, Opt. xvgr/xmgr file -px pullx.xvg Output, Opt. xvgr/xmgr file -pf pullf.xvg Output, Opt. xvgr/xmgr file -mtx nm.mtx Output, Opt. Hessian matrix -dn dipole.ndx Output, Opt. Index file Option Type Value Description -- -[no]h bool yes Print help info and quit -[no]version bool no Print version info and quit -niceint0 Set the nicelevel -deffnm string Set the default filename for all file options -xvg enum xmgrace xvg plot formatting: xmgrace, xmgr or none -[no]pd bool no Use particle decompostion -dd vector 0 0 0 Domain decomposition grid, 0 is optimize -npmeint-1 Number of separate nodes to be used for PME, -1 is guess -ddorder enum interleave DD node order: interleave, pp_pme or cartesian -[no]ddcheck bool yes Check for all bonded interactions with DD -rdd real 0 The maximum distance for bonded interactions with DD (nm), 0 is determine from initial coordinates -rconreal 0 Maximum distance for P-LINCS (nm), 0 is estimate -dlb enum autoDynamic load balancing (with DD): auto, no or yes -dds real 0.8 Minimum allowed dlb scaling of the DD cell size -gcomint-1 Global communication frequency -[no]v bool no Be loud and noisy -[no]compact bool yes Write a compact log file -[no]seppot bool no Write separate V and dVdl terms for each interaction type and node to the log file(s) -pforce real -1 Print all forces larger than this (kJ/mol nm) -[no]reprod bool no Try to avoid optimizations that affect binary reproducibility -cpt real
[gmx-users] pull code problem: between protofilaments
Dear Shilpi: Can you use something like this? pull = umbrella pull_geometry= position pull_dim = N N Y pull_vec1= 0 0 0 pull_start = no pull_ngroups = 1 pull_group0 = PRO-1 pull_pbcatom0= set this or use 0 pull_group1 = PRO-2 pull_pbcatom1= set this or use 0 pull_init1 = 0 0 set initial value here for z-axis pull_rate1 = 0 pull_k1 = 500.0 pull_nstxout = 500 pull_nstfout = 500 The above is for umbrella sampling. If you want to do continuous pulling, then: pull_start = yes pull_rate1 = set the rate ### Also: Next time you post, please provide more specifics. For example, I suggested a .mdp file in specifics to you above and I bet it would have been harder for you to guess what I meant if I had just told you the general idea instead of pasting some .mdp options. Likewise, your initial post would have been clearer it you had copied and pasted the .mdp pull code section that you tried to use. Chris. -- original message -- Dear Gmx users, I am studying the interaction between the tubulin protofilaments arranged in parallel. For this operation, I have considered a tetramer and a dimer from two protofilaments respectively (say PRO-1 and PRO-2), tetramer in PRO-1 and a dimer in PRO-2. I want to move the dimer of PRO-2 over the tetramer of PRO-1 along the length of protofilaments in one axis only, keeping the PRO-1 fixed to its original position. I tried by assuming tetramer as 'reference group' and the dimer as 'pull group' in pull code but the system crashed. I have succeeded in separating two dimers in Z-axis by using 'distance' geometry. But this case is quite different, as the pulling is not face-to-face but rather a sliding movement over another protofilament. Here, the COM distance between the pull group (dimer of PRO-2) and reference group (tetramer of PRO-1) first decreases and then increases while it moves. How can I simulate this operation by using pull code? Thanks, best regards, Shilpi Chaurasia -- 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] check point file corrupted - umbrella sampling
I am somewhat confused. When I do umbrella sampling, I use a bash script to set up a number of separate restrained simulations. Thus I would get one .cpt file for each simulation. What are you doing that in addition to these .cpt files you also get a state.cpt and state_perv.cpt? Perhaps you had better start at the beginning and explain what you did. Nevertheless, if your .cpt file is indeed corrupt, then you can not use it. You must extract a .gro from the final frame and use it (along with your .trr and .edr if you like) to generate a new .tpr file that will start a new run from where you left off. Chris. -- original message -- Hello Chris, Thank you very much for your reply. I have done following test as you suggested to check whether the checkpoint file are corrupted. 1) Yes, I have used same version of gromacs for mdrun and to produce .cpt file. 2) when i do gmxcheck -f 0.cpt # i get the same error. Fatal error: Start of file magic number mismatch, checkpoint file has 1993, should be 171817 The checkpoint file is corrupted or not a checkpoint file 3) When I do mdrun -s umbrella0.tpr # it runs without any error. 4) My mdruns has not outputted anything like 0_prev.cpt, i just have 0.cpt. I had run US for 24 windows, each for 10ns and produced 24 different cpt files. Along with these 24 different .cpt files, in the end it also produced state.cpt and state_prev.cpt. When i did gmxcheck -f state.cpt# i get Checking file state.cpt # Atoms 51895 Last frame -1 time 1.000 Item#frames Timestep (ps) Step 1 Time 1 Lambda 1 Coords 1 Velocities 1 Forces 0 Box 1 When i did gmxcheck -f state_prev.cpt # i get Checking file state_prev.cpt # Atoms 51895 Last frame -1 time 9730.730 Item#frames Timestep (ps) Step 1 Time 1 Lambda 1 Coords 1 Velocities 1 Forces 0 Box 1 I just have only one md.log file in the end, i haven't outputted md.log file for each window. Kind Regards, Chetan. From: gmx-users-boun...@gromacs.org [gmx-users-boun...@gromacs.org] On Behalf Of Chris Neale [chris.ne...@utoronto.ca] Sent: 19 September 2011 17:59 To: gmx-users@gromacs.org Subject: [gmx-users] check point file corrupted - umbrella sampling I have used tpbconv -until to extend a US simulation without error, and I assume that tpbconv -extend is also fine for US simulations. 1. Are you using the same version of gromacs mdrun as you used to produce your .cpt file initially? 2. what happens when you run gmxcheck on the .cpt file? 3. what happens when you run mdrun with the umbrella0.tpr file? 4. what happens when you use 0_prev.cpt? What I am getting at is that you have stated that the checkpoint file is corrupted in your subject line, but you have not provided any evidence of that. I do appreciate that you have pasted your exact commands. The next step is to see if there is any way to test your hypothesis that the checkpoint file is corrupted. Chris. -- original message -- Hi, I am using Gromacs 4.5.3 and want to extend my umbrella sampling runs. But i get an errror while running the mdrun step of extension command, saying Start of file magic number mismatch, checkpoint file has 1993, should be 171817. The checkpoint file is corrupted or not a checkpoint file Below are the commands I have used to run the umbrella sampling run grompp -f md_umbrella.mdp-c npt0.gro-t npt0.cpt-p topol.top-n index.ndx-o umbrella0.tpr mdrun-pf pullf-umbrella0.xvg-px pullx-umbrella0.xvg-s umbrella0.tpr-o 0.cpt-c 0.gro-x 0.xtc I used the below commands for extension: tpbconv -s umbrella0.tpr -extend 2 -o umbrella_new.tpr mdrun -s umbrella_new.tpr -cpi 0.cpt -pf pullf-umbrella_new.xvg -px pullx-umbrella_new.xvg -o new.cpt -c new.gro-x new.xtc (at this step i get the error) I have not touched the cpt files. Its the same as it was outputed. Please can I know what might have gone wrong. Kind Regards, chetan -- 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] check point file corrupted - umbrella sampling
I have used tpbconv -until to extend a US simulation without error, and I assume that tpbconv -extend is also fine for US simulations. 1. Are you using the same version of gromacs mdrun as you used to produce your .cpt file initially? 2. what happens when you run gmxcheck on the .cpt file? 3. what happens when you run mdrun with the umbrella0.tpr file? 4. what happens when you use 0_prev.cpt? What I am getting at is that you have stated that the checkpoint file is corrupted in your subject line, but you have not provided any evidence of that. I do appreciate that you have pasted your exact commands. The next step is to see if there is any way to test your hypothesis that the checkpoint file is corrupted. Chris. -- original message -- Hi, I am using Gromacs 4.5.3 and want to extend my umbrella sampling runs. But i get an errror while running the mdrun step of extension command, saying Start of file magic number mismatch, checkpoint file has 1993, should be 171817. The checkpoint file is corrupted or not a checkpoint file Below are the commands I have used to run the umbrella sampling run grompp -f md_umbrella.mdp-c npt0.gro-t npt0.cpt-p topol.top-n index.ndx-o umbrella0.tpr mdrun-pf pullf-umbrella0.xvg-px pullx-umbrella0.xvg-s umbrella0.tpr -o 0.cpt-c 0.gro-x 0.xtc I used the below commands for extension: tpbconv -s umbrella0.tpr -extend 2 -o umbrella_new.tpr mdrun -s umbrella_new.tpr -cpi 0.cpt -pf pullf-umbrella_new.xvg -px pullx-umbrella_new.xvg -o new.cpt -c new.gro-x new.xtc (at this step i get the error) I have not touched the cpt files. Its the same as it was outputed. Please can I know what might have gone wrong. Kind Regards, chetan -- 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] atom type OXT
You still haven't provided all of the relevant information. For
example: what ff are you using? But please don't just sent that
information now... better for you to read about how to develop a good
post and next time be suer to include enough information that we could
reproduce the problem if we tried.
In any event, the OXT atom is not recognized by your force field.
Depending on how many oxygen atoms you have in your final LEU residue,
you should rename it as you see based on either ff*.rtp or ff*-c.tdb.
I would personally keep only a single oxygen in your c-term LEU and
rename it O. You can add back the original coordinates by hand to
the .gro if you wish. I am not sure if pdb2gmx can handle being
supplied with the coordinates of the heavy atoms for the termini or
not. In any event, this will let pdb2gmx produce the .top/.itp for you.
I am a little surprised that you are running into this because
xlateat.dat should allow renaming of OXT -- O automatically but
perhaps the problem is that you are supplying both O and OXT? (but
then again you didn't provide the coordinates of the c-terminal
residue).
To be honest I have not struggled with this much because I have always
just done as I suggested to you above and it works fine. I guess it's
possible that you are using a very old version of gromacs (but then
again you didn't say what version you are using).
Chris.
Message: 4
Date: Tue, 06 Sep 2011 22:35:09 +1000
From: Mark Abraham Mark.Abraham at anu.edu.au
Subject: Re: [gmx-users] atom type OXT
To: Discussion list for GROMACS users gmx-users at gromacs.org
Message-ID: 4E66137D.9040103 at anu.edu.au
Content-Type: text/plain; charset=ISO-8859-1; format=flowed
On 6/09/2011 11:39 AM, Sweta Iyer wrote:
Hi all,
I have been trying to pdbgmx my protein to obtain the gro and top files
as follows:
pdb2gmx -f ${MOL}.pdb -o ${MOL}.gro -p ${MOL}.top -ter -ignh
However, I get an error message that states:
Atom OXT in residue SER 29 was not found in rtp entry SER with 8 atoms
while sorting atoms
We don't know what termini you are choosing (or want), so it's hard to
help. See also
http://www.gromacs.org/Documentation/How-tos/Multiple_Chains
Mark
Hi,
I chose NH3+ as the start terminus at the first residue which is leucine
and COO- as my end terminus which is also on a leucine residue, which is
when it says it cant recognize atom type OXT on my last leucine residue.
--
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] Constraints not working in pull code (sometimes, sometimes not)
Are you 100% sure that you have the correct index group for your large box? Did you use a .ndx file that you prepared as input to mdrun or did you rely on the creation of a standard group? If you relied on the creation of a standard group, can you please try again with the larger box and a .ndx file that you created (preferably the exact one that you used with g_dist) and supplied to mdrun? If that doesn't fix the problem, then I'd consider filing a redmine issue. Chris. -- original message -- Dear all, I have reported a week ago update about my troubles with octanol/water slab and pulling simulations, where COM code starts to report a wrong distances between centre-of-masses between two groups when bigger box (see bellow) is used. So far, smaller box cured our actual calculation problems, but I am still curious why g_dist and pullx.xvg are giving so different values, when I enlarge the simulation box? Best wishes Karel Berka Original mail down there Dear all, We have tested more strange behavior of constraints from pull code for calculation of free energy differences small molecule on DOPC membrane and octanol slab. As I have reported previously, while the errors in free energy differences on DOPC were rather small, the errors on octanol were strangely high. We have prepared smaller slabs (as Chris Neale suggested) of octanol with comparable size of box and we have also tried to analyze the position of centres either via pullx.xvg provided by pull code, but also by g_dist utility and to our surprise, size matters. A lot. (Same groups were used in mdp for pull and for g_dist analysis, only distances in z-axis is shown) DOPC, 128 molecules, box size 6.2 6.4 8.3 , position of small molecule on the edge of membrane Time(ps) pullf.xvgdist.xvg 0.00-0.12416-0.12436 1.00-0.12416-0.12433 2.00-0.12416-0.12428 3.00-0.12416-0.12430 4.00-0.12416-0.12430 5.00-0.12416-0.12432 6.00-0.12416-0.12434 7.00-0.12416-0.12449 As can be seen, the variation of distances reported by these two programs is not perfect, but it is rather ok with precision of about 0.001 nm. Octanol, 957 molecules, box size 5.7 5.9 13.7, position is similar on the edge of slab, errors in free energy profiles are mainly here. Time(ps)pullf.xvgdist.xvg 0.00-2.81855-3.41133 1.00-2.81076-3.44213 2.00-2.82005-3.27949 3.00-2.81856-3.35097 4.00-2.82016-3.50378 5.00-2.81849-3.63261 6.00-2.81855-3.60058 7.00-2.81870-3.80251 8.00-2.81859-3.32790 9.00-2.81849-3.24329 10.00-2.81855-3.28394 21.00-2.81855-4.02524 Here, the distances reported by pull code and by g_dist are completely different with differences in about 1 nm! Visual inspection in VMD says that small molecule is highly mobile a actual distances are reported rightly by g_dist program. Strangely enough, when we put distance measured by g_dist as a start into the pull code without guessing the distances, but this also did not worked, since the molecule tried to be at some completely different position and the simulation crashed as the molecule speeded up. Octanol, 475 molecules, box size 5.5 5.5 9.5, position also on the edge of the slab Time(ps) pullf.xvgdist.xvg 0.00-0.49996-0.49989 1.00-0.49996-0.49988 2.00-0.49996-0.49981 3.00-0.49996-0.49981 4.00-0.49996-0.50001 5.00-0.49996-0.49987 6.00-0.49996-0.49986 7.00-0.49996-0.49994 8.00-0.49997-0.49988 9.00-0.49996-0.49980 10.00-0.49996-0.49988 Smaller box cured miraculously the problem, since here the difference is similarly small as in DOPC case. However just to remind another ache we had with octanol slab, in case of octanol, the size of box in xy plane had to be held constant as the box narrowed thorough the whole simulation when in anisotropic NPT ensemble. Why g_dist and pullx show so different results in large box, I still do not understand. Any suggestions? And possibly any suggestions whether our box narrowing in case of octanol is common or do we have another devil hidden there? Best wishes -- Zdraví skoro zdravý Karel Krápník Berka -- 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] will united-atom POPC forcefield mater in pH 4
Zhijun: This is a question for the gmx-users list, so I will address it there. Please keep this type of question on the users list. 1.Are hydrogen ion H+ probably bind to POPC's phosphate O atom when I randomly assign it? That sounds a lot like something that you can test. Na+ ions do not form long-lived interactions like this during simulation (spending much more time in bulk water than in particular interactions with individual lipids), so I doubt it. 2.Will hydrogen ion (so small) channel the POPC lipid, while a full-atom is not? Your question is not well formed. Do you mean will the hydrogen ion pass between the lipids, across the bilayer in united-atom but not all-atom simulations? Again, that sounds like something that you can test. It would eventually cross both types of bilayers (in the limit of infinite simulations) but no, I doubt that you would observe this in obtainable simulation timescales. Chris. -- original message -- Hi,dear NAMDers I am a gromacs user but recently run a simultion at pH 4, which might protonate the protein embedded in POPC(a phosphatidylcholine) membrane.But I'm using a united-atom model to POPC develop by Erik Lindahl and Peter Kasson(2010), base on Berger lipids ported to Gromacs/AMBER forcefields(1997). I roughly calculates what's pH 4 means:10^4 mol/dm^3=1/18000nm^3, about one H+ per 27nm cubic box. There are two things I am not quiet sure : 1.Are hydrogen ion H+ probably bind to POPC's phosphate O atom when I randomly assign it? 2.Will hydrogen ion (so small) channel the POPC lipid, while a full-atom is not? Do you have some advise for how to do? Or a full-atom model for POPC better? -- 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] Half double pair list method in GROMACS
I think it's reasonable for you to post this on-list. It's not your method and you got help here and somebody else might want the same help in the future. Chris. -- original message -- ABEL Stephane 175950 Stephane.ABEL at cea.fr Thu Sep 1 22:22:27 CEST 2011 * Previous message: [gmx-users] pdb2gmx response time * Next message: [gmx-users] order parameter * Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] Chris, Thank you for your confirmation. I did the changes. I am currently doing some tests, I will send you a feedback about the results off-list (if you want) shortly. A bientôt Stéphane Message: 1 Date: Thu, 01 Sep 2011 14:38:53 -0400 From: chris.neale at utoronto.ca Subject: [gmx-users] Half double pair list method in GROMACS To: gmx-users at gromacs.org Message-ID: 20110901143853.537ln0dj94w4wc4c at webmail.utoronto.ca Content-Type: text/plain; charset=ISO-8859-1; DelSp=Yes; format=flowed It's a typo, but it's in the discussion and not in the do this part of the method so I decided not to mention it. I don't see another question in this post, so I hope that you have figured things out. Note that I have never tested the exact implementation that I suggested in that April email. It seems like it should work just fine, but it is numerically different than the OPLS/Berger combination so there is no way to be sure until you check the energies as I suggested. I'd be interested to have you report back on the energy matching once you have done the tests. Good luck, Chris. -- original message -- Hi Chris, Sorry to repost the same question, but I have really tested your method the last few weeks. My question about the gen-pairs directive come from the fact that I have read a message from you http://lists.gromacs.org/pipermail/gmx-users/2006-September/023761.html Where you detailed how to use the Berger and OPLS force fields together in the same system. By reading carefully the meaning of the gen-pairs directive, I found several errors in my force field. BYW in your previous message http://lists.gromacs.org/pipermail/gmx-users/2011-April/060839.html In the step 3, I think there is a typo, it is values/10 instead of values/12. Am I right? Thank you again Stephane -- 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] GROMACS 4.5.4 keep crashing all the time
I am glad that the pressure coupling intervals have been identified as a source of instability for poorly equilibrated systems as I was unaware of that. Still, the fact that the SD integrator also solves the problem also suggests that this is simply a poorly equilibrated system. I am not sure why PME would run fine and reaction field would give lincs warnings, but then again I have no experience with using a reaction field. Chris. On 1/09/2011 6:32 PM, Itamar Kass wrote: Hi Tsjerk, Thanks for the help, it actually worked. When nstpcouple is set to 1m the system can be equilibrated (NPT) without LINCS error. I hadn't thought about it as I never use NVT (unless doing FE calculations). Equilibrating with NVT before NPT can be wise when the system starts far from equilibrium. Is this mean the 4.5.4 is more sensitive the 4.0.7? Is this effect the data collected till now? If this is the case, why 5 ns simulations done with 4.0.7 crashed when extended it using 4.5.4? IIRC 4.0.x and 4.5.x have different mechanisms for deciding how often to do global communication for things like T and P coupling. Mostly you can get away with the same kind of approximation one uses with twin-range neighbour lists, periodic neighbour list updates, RESPA, etc. where slowly-varying quantities don't have to be recalculated every step. However during equilibration (and that includes the transition from 4.0.x to 4.5.x) those assumptions need not be valid. So tuning the update frequency to be high during transitions is a good idea. Then relax them and see how you go. Mark Also, is this mean I can do my productive run using 4.5.4 with the default value of nstpcouple, it seems that setting it to 1 greatly increases the computational time. To the best of my knowledge, in prior version nstpcouple was set to 1 by default. Cheers, Itamar On 01/09/2011, at 4:47 PM, Tsjerk Wassenaar wrote: Hi Itamar, I haven't really followed the discussion and I'm a bit too lazy to look it all up now ;) But have you tried setting the nst parameters to 1 (except for output). Especially nstpcouple. Note that nstpcouple=1 requires nstlist=1 and nstcalcenergy=1. If that solves the problem, you may need to extend your equilibration a bit, first relaxing NVT, followed by NPT with nstpcouple=1, thereafter equilibrating using productions conditions. It it solves it, maybe the option should be renamed nstptrouble :p Hope it helps, Tsjerk -- 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] Non-zero total charge
6.30e+01 = 63 -- original message -- My system had a no zero total charge: System has non-zero total charge: 6.30e+01 I used genion to neutralize the system by adding 6 CL ions. After updating the topology file, the system still seems to have the same problem. It still has a non-zero total charge: System has non-zero total charge: 5.70e+01 Please let me how i can rectify the situation. Many thanks, Munishika -- next part -- An HTML attachment was scrubbed... URL: http://lists.gromacs.org/pipermail/gmx-users/attachments/20110901/58ff645a/attachment.html -- 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] GROMACS 4.5.4 keep crashing all the time
That all makes sense Tsjerk. I wonder if mdrun terminations based on LINCS warnings should come with an additional message to explain that one may try running for a while with nstpcouple=1. Also, I'm still a little curious about a question that Itamar asked a few posts ago: If this is the case, why 5 ns simulations done with 4.0.7 crashed when extended it using 4.5.4? Mark provided a good explanation about how this could be possible, but I've never seen lincs warnings or systems blowing up after 5 ns of equilibration. I fully realize that it may take even us of simulation to properly equilibrate statistical properties, but in my experience it is far outside of ordinary to require 5 ns of equilibration to avoid systems blowing up. Chris. -- original message -- Hi Chris, I can imagine that the pressure scaling has a more profound effect on the 'visible' surroundings if a cut-off is used, while this will not be the case when using PME. So the shock for an atom when the coordinates are adjusted can be expected to be greater with cut-off based methods, rendering such simulations less stable. As for SD, probably that causes sufficient damping of jitter introduced due to pressure coupling for it not to propagate and cause problems. But those are just my two cents (about 2.8 dollar cents with current rates :p). Cheers, Tsjerk On Thu, Sep 1, 2011 at 2:41 PM, chris.neale at utoronto.ca wrote: I am glad that the pressure coupling intervals have been identified as a source of instability for poorly equilibrated systems as I was unaware of that. Still, the fact that the SD integrator also solves the problem also suggests that this is simply a poorly equilibrated system. I am not sure why PME would run fine and reaction field would give lincs warnings, but then again I have no experience with using a reaction field. Chris. On 1/09/2011 6:32 PM, Itamar Kass wrote: Hi Tsjerk, Thanks for the help, it actually worked. When nstpcouple is set to 1m the system can be equilibrated (NPT) without LINCS error. I hadn't thought about it as I never use NVT (unless doing FE calculations). Equilibrating with NVT before NPT can be wise when the system starts far from equilibrium. Is this mean the 4.5.4 is more sensitive the 4.0.7? Is this effect the data collected till now? If this is the case, why 5 ns simulations done with 4.0.7 crashed when extended it using 4.5.4? IIRC 4.0.x and 4.5.x have different mechanisms for deciding how often to do global communication for things like T and P coupling. Mostly you can get away with the same kind of approximation one uses with twin-range neighbour lists, periodic neighbour list updates, RESPA, etc. where slowly-varying quantities don't have to be recalculated every step. However during equilibration (and that includes the transition from 4.0.x to 4.5.x) those assumptions need not be valid. So tuning the update frequency to be high during transitions is a good idea. Then relax them and see how you go. Mark Also, is this mean I can do my productive run using 4.5.4 with the default value of nstpcouple, it seems that setting it to 1 greatly increases the computational time. To the best of my knowledge, in prior version nstpcouple was set to 1 by default. Cheers, Itamar On 01/09/2011, at 4:47 PM, Tsjerk Wassenaar wrote: Hi Itamar, I haven't really followed the discussion and I'm a bit too lazy to look it all up now ;) But have you tried setting the nst parameters to 1 (except for output). Especially nstpcouple. Note that nstpcouple=1 requires nstlist=1 and nstcalcenergy=1. If that solves the problem, you may need to extend your equilibration a bit, first relaxing NVT, followed by NPT with nstpcouple=1, thereafter equilibrating using productions conditions. It it solves it, maybe the option should be renamed nstptrouble :p Hope it helps, Tsjerk -- gmx-users mailing listgmx-users at 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 thewww interface or send it to gmx-users-request aGROMACS 4.5.4 keep crashing all the time Tsjerk Wassenaar tsjerkw at gmail.com Thu Sep 1 15:48:40 CEST 2011 * Previous message: [gmx-users] GROMACS 4.5.4 keep crashing all the time * Next message: [gmx-users] Non-zero total charge * Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] Hi Chris, I can imagine that the pressure scaling has a more profound effect on the 'visible' surroundings if a cut-off is used, while this will not be the case when using PME. So the shock for an atom when the coordinates are adjusted can be expected to be greater with cut-off based methods, rendering such simulations less stable. As for SD, probably that causes sufficient damping of jitter introduced due to pressure coupling for it not to propagate and cause problems. But
[gmx-users] half double pair list method in GROMACS
Dear Stephane: We discussed this in April: http://lists.gromacs.org/pipermail/gmx-users/2011-April/060839.html At that time I also provided a method for you to verify your files (and the method in general). It is possible for you to answer your gen-pairs question by looking into the manual and reading about pairs, gen-pairs, and pairtypes. I think that this is one case where you would benefit more from fully understanding how these parts work than from a direct answer to your question. If you are having problems, please provide a whole bunch more information on the problems that you are seeing. If, on the other hand, you are just looking for somebody other than me to comment on the accuracy of the April post, then that is perfectly fine with me, but you should state that. Chris. -- original message -- Dear All, I try to apply the half double pair list method for a system containing a glycolipid surfactant and a peptide modeled with the GLYCAM and AMBER99SB force field. Briefly what I did : 1- I have changed the forcefield.itp like this [ defaults ] ; gen_pairs set explicitly --- gen-pairs = no ; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ ;1 2 yes 0.5 0.8333 --- for AMBER99SB only 1 2no 1.0 0.16 --- for both GLYCAM et AMBER99SB ;1 2yes 1.0 1.0 -- for GLYCAM only #include ffnonbonded_mod.itp #include ffbonded.itp 2- For the surfactant, I have calculated the pair-types interactions manually with the comb-rule = 2 and divided the values /6 and added these values in [pairtypes] section in the ffnonbonded_mod.itp files 3- For the peptide, I have calculated the pair-types interactions manually comb-rule = 2 and divided the values /10 and added these values in [pairtypes] section in the ffnonbonded_mod.itp files. 4- In the surfactant topology file I have repeated 6 times the [pairs] directives 0.16*6 ~=10 5 - In the peptide topology file I have repeated 5 times the [pairs] directives 0.16*5 ~= 0.8333 Is it correct ? However I have a little doubt about gen-pairs directive should I have set it to no or yes. in a previous message with a similar problem, the gen directive was set yes http://lists.gromacs.org/pipermail/gmx-users/2006-September/023761.html Thank you for your help Stephane -- 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] order parameter
In addition to Justin's comment about g_order being incorrect for unsaturated carbons, you won't ever get numerical agreement with g_order, even for saturated carbons, because g_order does not use your explicit hydrogens. g_order uses the positions of the carbons to rebuild the hydrogen positions assuming perfect tetrahedral geometry, which you certainly don't have in all cases when you use real hydrogens. There is a VMD tcl script that you can download that will compute order parameters from real hydrogen atoms. perhaps you should try to compare to that. Chris. -- original message -- Hi, Thanks lot for replying. I am doing all-atom simulation. I am doing the PBC before finding the angle. also I normalise the vectors before finding the angle. Yes I have checked that the formula I am using for the order parameter is correct. I am doing the averaging correctly. 1. I take the carbons in each tail ( I neglect the 1st and the last as GROMACS does) , then I find the Hydrogens associated with it. 2. Then I do the PBC , normalise them and then take the angle, then calculate the order parameter. 3. finally I average them over the frames. I have gone through the procedure and still I am not getting the same profile as GROMACS gives. Is there anything else that I need to include in my calculations? Ramya From: gmx-users-bounces at gromacs.org [gmx-users-bounces at gromacs.org] on behalf of chris.neale at utoronto.ca [chris.neale at utoronto.ca] Sent: Wednesday, August 31, 2011 8:00 PM To: gmx-users at gromacs.org Subject: [gmx-users] order parameter Dear Ramya: Are you simulating all-atom lipids (with explicit hydrogen atoms on the acyl chain)? If not, then you missed a step in your description of what you have done (g_order, for example, ignores explicit hydrogen atoms so that it can act on united atom lipids). Not sure why PBC would be your step #3, after your step #2 was to find the angle. I suggest that you simply run trjconv -pbc mol on your trajectory file before you process it and then you no longer need to worry about PBC in your custom analysis tool. Once you have the angle, you must average it correctly. The equation is available in most papers that describe order parameters and is listed as a comment at the top of the gmx_order.c source file (in version 4.0.7 at least). If you want to get more help on your procedure after you have worked on this for a while, I suggest laying out your procedure very specifically. Your previous post, for example, was pretty loose with terminology when you described your method and there is quite a bit that one must assume. Chris. -- original message -- Hi, I am trying to write a code for Deuterium order parameter of DOPC lipid. I went through the code in gmx_order.c, I did the following, 1. I took the carbons in the chain, and found its neighbors. 2. Took the bilayer normal and found the angle between the bilayer normal and the ?CH molecular axis. 3. Took care of the periodic boundary conditions since I use NPT ensemble. But the code in gmx_order.c in GROMACS tries to do a lot of things other than this, as I don?t know C or C++ language that it is using, I don?t know what else I am supposed to include. Can someone please help me? Ramya -- -- 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] Half double pair list method in GROMACS
It's a typo, but it's in the discussion and not in the do this part of the method so I decided not to mention it. I don't see another question in this post, so I hope that you have figured things out. Note that I have never tested the exact implementation that I suggested in that April email. It seems like it should work just fine, but it is numerically different than the OPLS/Berger combination so there is no way to be sure until you check the energies as I suggested. I'd be interested to have you report back on the energy matching once you have done the tests. Good luck, Chris. -- original message -- Hi Chris, Sorry to repost the same question, but I have really tested your method the last few weeks. My question about the gen-pairs directive come from the fact that I have read a message from you http://lists.gromacs.org/pipermail/gmx-users/2006-September/023761.html Where you detailed how to use the Berger and OPLS force fields together in the same system. By reading carefully the meaning of the gen-pairs directive, I found several errors in my force field. BYW in your previous message http://lists.gromacs.org/pipermail/gmx-users/2011-April/060839.html In the step 3, I think there is a typo, it is values/10 instead of values/12. Am I right? Thank you again Stephane Dear Stephane: We discussed this in April: http://lists.gromacs.org/pipermail/gmx-users/2011-April/060839.html At that time I also provided a method for you to verify your files (and the method in general). It is possible for you to answer your gen-pairs question by looking into the manual and reading about pairs, gen-pairs, and pairtypes. I think that this is one case where you would benefit more from fully understanding how these parts work than from a direct answer to your question. If you are having problems, please provide a whole bunch more information on the problems that you are seeing. If, on the other hand, you are just looking for somebody other than me to comment on the accuracy of the April post, then that is perfectly fine with me, but you should state that. Chris. -- original message -- Dear All, I try to apply the half double pair list method for a system containing a glycolipid surfactant and a peptide modeled with the GLYCAM and AMBER99SB force field. Briefly what I did : 1- I have changed the forcefield.itp like this [ defaults ] ; gen_pairs set explicitly --- gen-pairs = no ; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ ;1 2 yes 0.5 0.8333 --- for AMBER99SB only 1 2no 1.0 0.16 --- for both GLYCAM et AMBER99SB ;1 2yes 1.0 1.0 -- for GLYCAM only #include ffnonbonded_mod.itp #include ffbonded.itp 2- For the surfactant, I have calculated the pair-types interactions manually with the comb-rule = 2 and divided the values /6 and added these values in [pairtypes] section in the ffnonbonded_mod.itp files 3- For the peptide, I have calculated the pair-types interactions manually comb-rule = 2 and divided the values /10 and added these values in [pairtypes] section in the ffnonbonded_mod.itp files. 4- In the surfactant topology file I have repeated 6 times the [pairs] directives 0.16*6 ~=10 5 - In the peptide topology file I have repeated 5 times the [pairs] directives 0.16*5 ~= 0.8333 Is it correct ? However I have a little doubt about gen-pairs directive should I have set it to no or yes. in a previous message with a similar problem, the gen directive was set yes http://lists.gromacs.org/pipermail/gmx-users/2006-September/023761.html Thank you for your help Stephane -- 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] GROMACS 4.5.4 keep crashing all the time
That seems possible Mark. I had actually assumed that Itamar extracted a .gro from the 4.0.7 simulations and created a new .tpr with 4.5.4 grompp. If that was the case, then I would still be surprised that it had instability problems. One sees a lot of charmm papers where people do slow heating, but my experience with gromacs is that not too many people worry about maintaining the velocities from their NVT simulations for their NPT simulations (or restrained - unrestrained) and that it works out just fine. On the other hand, I don't see too many gromacs reaction field studies either, so perhaps it all makes sense to those who use these methods. Thanks again, Chris. On 2/09/2011 12:59 AM, chris.neale at utoronto.ca wrote: That all makes sense Tsjerk. I wonder if mdrun terminations based on LINCS warnings should come with an additional message to explain that one may try running for a while with nstpcouple=1. That tip might be a good thing for the wiki page on that. Also, I'm still a little curious about a question that Itamar asked a few posts ago: If this is the case, why 5 ns simulations done with 4.0.7 crashed when extended it using 4.5.4? Mark provided a good explanation about how this could be possible, but I've never seen lincs warnings or systems blowing up after 5 ns of equilibration. I fully realize that it may take even us of simulation to properly equilibrate statistical properties, but in my experience it is far outside of ordinary to require 5 ns of equilibration to avoid systems blowing up. It's hard to say without more detail of how the extension occurred, and knowing how much ensemble data got lost. It's still conceivable some horrible mismatch occurred (e.g. virial over-written by some other data), but not really worth exploring properly. I'd just expect to have to re-equilibrate upon changing code version, and just not attempt such an extension. Mark Chris. -- original message -- Hi Chris, I can imagine that the pressure scaling has a more profound effect on the 'visible' surroundings if a cut-off is used, while this will not be the case when using PME. So the shock for an atom when the coordinates are adjusted can be expected to be greater with cut-off based methods, rendering such simulations less stable. As for SD, probably that causes sufficient damping of jitter introduced due to pressure coupling for it not to propagate and cause problems. But those are just my two cents (about 2.8 dollar cents with current rates :p). Cheers, Tsjerk On Thu, Sep 1, 2011 at 2:41 PM, chris.neale at utoronto.ca wrote: I am glad that the pressure coupling intervals have been identified as a source of instability for poorly equilibrated systems as I was unaware of that. Still, the fact that the SD integrator also solves the problem also suggests that this is simply a poorly equilibrated system. I am not sure why PME would run fine and reaction field would give lincs warnings, but then again I have no experience with using a reaction field. Chris. On 1/09/2011 6:32 PM, Itamar Kass wrote: Hi Tsjerk, Thanks for the help, it actually worked. When nstpcouple is set to 1m the system can be equilibrated (NPT) without LINCS error. I hadn't thought about it as I never use NVT (unless doing FE calculations). Equilibrating with NVT before NPT can be wise when the system starts far from equilibrium. Is this mean the 4.5.4 is more sensitive the 4.0.7? Is this effect the data collected till now? If this is the case, why 5 ns simulations done with 4.0.7 crashed when extended it using 4.5.4? IIRC 4.0.x and 4.5.x have different mechanisms for deciding how often to do global communication for things like T and P coupling. Mostly you can get away with the same kind of approximation one uses with twin-range neighbour lists, periodic neighbour list updates, RESPA, etc. where slowly-varying quantities don't have to be recalculated every step. However during equilibration (and that includes the transition from 4.0.x to 4.5.x) those assumptions need not be valid. So tuning the update frequency to be high during transitions is a good idea. Then relax them and see how you go. Mark Also, is this mean I can do my productive run using 4.5.4 with the default value of nstpcouple, it seems that setting it to 1 greatly increases the computational time. To the best of my knowledge, in prior version nstpcouple was set to 1 by default. Cheers, Itamar On 01/09/2011, at 4:47 PM, Tsjerk Wassenaar wrote: Hi Itamar, I haven't really followed the discussion and I'm a bit too lazy to look it all up now ;) But have you tried setting the nst parameters to 1 (except for output). Especially nstpcouple. Note that nstpcouple=1 requires nstlist=1 and nstcalcenergy=1. If that solves the problem, you may need to extend your equilibration a bit, first relaxing NVT, followed by NPT with nstpcouple=1, thereafter equilibrating using productions conditions. It it
[gmx-users] order parameter
Dear Ramya: Are you simulating all-atom lipids (with explicit hydrogen atoms on the acyl chain)? If not, then you missed a step in your description of what you have done (g_order, for example, ignores explicit hydrogen atoms so that it can act on united atom lipids). Not sure why PBC would be your step #3, after your step #2 was to find the angle. I suggest that you simply run trjconv -pbc mol on your trajectory file before you process it and then you no longer need to worry about PBC in your custom analysis tool. Once you have the angle, you must average it correctly. The equation is available in most papers that describe order parameters and is listed as a comment at the top of the gmx_order.c source file (in version 4.0.7 at least). If you want to get more help on your procedure after you have worked on this for a while, I suggest laying out your procedure very specifically. Your previous post, for example, was pretty loose with terminology when you described your method and there is quite a bit that one must assume. Chris. -- original message -- Hi, I am trying to write a code for Deuterium order parameter of DOPC lipid. I went through the code in gmx_order.c, I did the following, 1. I took the carbons in the chain, and found its neighbors. 2. Took the bilayer normal and found the angle between the bilayer normal and the ?CH molecular axis. 3. Took care of the periodic boundary conditions since I use NPT ensemble. But the code in gmx_order.c in GROMACS tries to do a lot of things other than this, as I don?t know C or C++ language that it is using, I don?t know what else I am supposed to include. Can someone please help me? Ramya -- 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] GROMACS 4.5.4 keep crashing all the time
Itamar: We really are trying to help. I think that perhaps you don't grasp how difficult it is to help without being able to access the simulation directly. Therefore we have ideas and we ask you to do specific things that are going to move us toward a solution, either by finding answers or by ruling out possibilities. It is actually useful information to know that Sometimes it is the peptide N and H, like in the case of 981 and 982, and sometimes others ... but when it seems like you don't want to provide the requested information my first inclination is to give up on trying to help. At this point, there are a few unanswered old questions and I have some new questions. 1. Can you reproduce this with a water box? 2. Can you reproduce this with your protein in vacuum? 3. If neither 2 or 3, then can you step slowly from one of these systems toward your final system and identify the point at which the lincs warnings arise? 4. Do you get the warnings without Ca also, or just with Ca? 5. Can you reproduce this with the SD integrator? If you are really against trying this, then at least can you reproduce this with a single Berendsen temperature coupling group? 6. Can you reproduce this without using the reaction field? Either with PME or a simple cutoff? 7. Can you trace down the version of gromacs (between 4.0.7 and 4.5.4) where you start to see this warning? Chris. -- original message -- Hi Justin, I did repeat it using gen_val and running temperature the same, with no effect, it is still crash. I didn't replied point #6 because the atoms which triggers the LINCS are different between each try. Sometimes it is the peptide N and H, like in the case of 981 and 982, and sometimes others. In addition, there is no visible difference in dynamics between 4.0.7 and 4.5.4 I can find. Itamar On 01/09/2011, at 11:14 AM, Justin A. Lemkul wrote: Itamar Kass wrote: Hi Mark, I didn't had the time to do the SD yet, but serial run end with the same results. I didn't try water only system, as this is of no interest to me, but I will simplify the system later on. Being of interest to you and being a useful diagnostic may be different. It's important to rule out different variables to arrive at a solution, which I suspect is of interest to you. You also haven't addressed points #1 and #6 in Chris' message. -Justin Cheers, Itamar On 01/09/2011, at 10:51 AM, Mark Abraham wrote: On 1/09/2011 10:20 AM, Itamar Kass wrote: Hi Chris, Thanks for the email, I am sorry it took me some time to replay. I tried 4.5.4 again, now starting from a 5 ns simulations run using 4.0.7, and again the simulations had stopped after 1000 LINCS error (I can extend the simulations using 4.0.7). I know that gromacs stopped after 1000 LINCS, but this is usually a sign that something bad is going on in the system. OK. Chris suggested a number of other strategies that will help determine which aspect of 4.5.4 is behaving differently. How did those strategies work out? Mark Cheers, Itamar On 18/08/2011, at 12:03 PM, chris.neale at utoronto.ca wrote: OK, here's my last few ideas: 1. Please try to repeat this with gen_vel set to the same value as your temperature coupling 2. Can you reproduce this in serial? 3. Can you reproduce this with the sd integrator? 4. Can you reproduce this with a simpler system? protein in vacuum or just water or remove the ions, etc? 5. Take the output .gro from 4.0.7 that ran fine for X ns and run it under 4.5.4. Do you get the same lincs warnings? 6. Also, note that you are getting warnings and the run does not actually crash but just stops after too many warnings. So what are atoms 981 and 982? Does their motion look different in an important ways between the 4.0.7 and 4.5.4 trajectories? Chris. -- original message -- Hi Chris, thanks for the advice, I have to say I tried this as well without any success. Itamar On 18/08/2011, at 11:11 AM, chris.neale at utoronto.ca wrote: run an EM with flexible water. I often find that this is the only way to get a stable system. 500 steps of steep with define=-DFLEXIBLE (or different depending on your water model I think) should be enough. Chris. Hi Chris and Justin, On 18/08/2011, at 9:36 AM, Justin A. Lemkul wrote: -- gmx-users mailing listgmx-users at 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-request at gromacs.org. Can't post? Read http://www.gromacs.org/Support/Mailing_Lists - In theory, there is no difference between theory and practice. But, in practice, there is. - Jan L.A. van de Snepscheut === | Itamar Kass, Ph.D. | Postdoctoral Research Fellow | | Department of
[gmx-users] g_wham issue
Please provide very much more information. Unless somebody has run into the exact thing that you are describing it is currently impossible for us to help you. One thing I can suggest now is to construct histograms of the sampling from the coord.xvg files yourself with some scripting and plot them all together to see if you have some regions with no overlap (but even still, please address my first point). -- original message -- Hi all: I have been trying to figure out a problem in our umbrella sampling simulation. We are simulating a coarse grained protein passing through a coarse grained lipid bilayer. In the pull simulation, the protein was made to move from +z=6 nm of the box through the lipid layer at z=0 to z= -4 nm. For umbrella sampling 74 windows where chosen with a 0.2 nm. When we look at the histogram from the g_wham, we get good overlap before the protein touches the bilayer and after it has exited the bilayer. We also get the PMF curves for these regions. However, when the the protein is in contact, and barely touching the bilayer, the g_wham analysis gives 0.000E+0 for PMF and no histograms in the output. We have tried decreasing the window size from 0.2 to 0.1 nm but it does not help. I looked at the gmx-user list to find a possible solution to this issue but could not find any helpful clue. Could anyone point me to what should be done and what is possibly going wrong? Thanks, SNC -- next part -- An HTML attachment was scrubbed... URL: http://lists.gromacs.org/pipermail/gmx-users/attachments/20110825/03e56b1f/attachment-0001.html -- 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] Free Energy Integrator Selection
Try running for longer and look at the convergence of the temperatures over time. I suspect that what you are reporting is actually just statistical noise. I am not sure about using the sd integrator and defining multiple temperature coupling groups... If you define the same temperature for each then I suspect this is the exact same as defining the system as a single temperature coupling group (for the sd integrator). Since you can not extract dynamic information from your FE simulations anyway, the sd integrator is the best choice. Chris. -- original message -- Hi there, I am also doing FEP calculations and I am also using sd as integrator. The problem that I am having with this integrator is the temperature coupling (NTP). I am setting the targeted temperature of the system on 300K, for 2 Groups, Protein_and_ligand and Solvent. g_energy gives me 303K 301K respectably. If I use the md and noose-hover thermostat the results are 300K for both groups. For that reason I was thinking on changing to md. But afther reading the post http://lists.gromacs.org/pipermail/gmx-users/2007-February/025815.html I think that I shuold stick on sd. The questions are If by trying sd with target temperature 298K and get 300K as a result, is the simulation valid? Fabian, If you use sd the temperature of your system is the wanted one? Thanks in advance Marcelino Justin A. Lemkul jalemkul at vt.edu wrote: -- 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: parallel job crash for large system
Your density seems to be about 70% of what I would expect. Are you sure that this is not just a normal case of a poorly equilibrated system crashing? That matches with what you say about the density growing (although perhaps it has more to do with poor equilibration than with mixing, as you suggest)? In any event, I'd suggest simplifying your system and making it smaller to see of you can reproduce the problem with a system that will run quickly in serial. Chris. On 23/08/2011 8:44 AM, Dr. Vitaly V. Chaban wrote: In the below issue, the barostat is setup semiisotropically and works only along the long direction. The density of the system slowly grows due to mixing. If this can be useful Does a different barostat work? Mark On Mon, Aug 22, 2011 at 5:32 PM, Dr. Vitaly V. Chaban vvchaban at gmail.com wrote: We are running the system consisting of 84000 atoms in parallelepipedic box, 6x6x33nm. The starting geometry, etc are OK and evolution of trajectory is reasonable but after several hundred thousands of steps it suddenly crashes. Mysteriously, each time it crashes at different time-steps, but it always occurs. The parts of this system were equilibrated separately and did not crash. The system is not in equilibrium but without external forces. The Parrinello-Rahman barostat is turned on. The md.log does not show any problems, the PDB configurations are not written down before crash, the constaints are absent, the time-step is 1fs that is OK for separate components (in separate boxes). With serial gromacs, the error is not yet observed, but given the size the run is very slow. What can it be? Can it be somehow connected with the very (oblongated) box? Stdout below: 5000 steps, 5.0 ps. [exciton04:10256] *** Process received signal *** [exciton04:10256] Signal: Segmentation fault (11) [exciton04:10256] Signal code: Address not mapped (1) [exciton04:10256] Failing at address: 0x6c0ebf10 [exciton04:10257] *** Process received signal *** [exciton04:10257] Signal: Segmentation fault (11) [exciton04:10257] Signal code: Address not mapped (1) [exciton04:10257] Failing at address: 0x6378320 [exciton04:10253] *** Process received signal *** [exciton04:10253] Signal: Segmentation fault (11) [exciton04:10253] Signal code: Address not mapped (1) [exciton04:10253] Failing at address: 0x1bfbe110 [exciton04:10253] [ 0] /lib64/libpthread.so.0 [0x3402a0eb10] [exciton04:10253] [ 1] mdrun [0x66bb4d] [exciton04:10253] *** End of error message *** [exciton04:10255] *** Process received signal *** [exciton04:10255] Signal: Segmentation fault (11) [exciton04:10255] Signal code: Address not mapped (1) [exciton04:10255] Failing at address: 0x13dd139b0 [exciton04:10255] [ 0] /lib64/libpthread.so.0 [0x3402a0eb10] [exciton04:10255] [ 1] mdrun [0x66bb5e] [exciton04:10255] *** End of error message *** [exciton04:10256] [ 0] /lib64/libpthread.so.0 [0x3402a0eb10] [exciton04:10256] [ 1] mdrun [0x66bb6f] [exciton04:10256] *** End of error message *** [exciton04:10254] *** Process received signal *** [exciton04:10254] Signal: Segmentation fault (11) [exciton04:10254] Signal code: Address not mapped (1) [exciton04:10254] Failing at address: 0x13d2103b0 [exciton04:10254] [ 0] /lib64/libpthread.so.0 [0x3402a0eb10] [exciton04:10254] [ 1] mdrun [0x66bb5e] [exciton04:10254] *** End of error message *** [exciton04:10257] [ 0] /lib64/libpthread.so.0 [0x3402a0eb10] [exciton04:10257] [ 1] mdrun [0x66bb5e] [exciton04:10257] *** End of error message *** -- mpirun noticed that process rank 0 with PID 10253 on node exciton04 exited on signal 11 (Segmentation fault). -- 5 total processes killed (some possibly by mpirun during cleanup) The version is 4.0.7 used with OpenMPI. -- 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] problem with g_density -center
g_density is not compatible with constant pressure simulations. You
must modify it to construct the bins outward from the bilayer center
when doing NPT:
http://lists.gromacs.org/pipermail/gmx-users/2010-November/055651.html
Further, trjconv -center is misleading. I actually lost a lot of time
thinking that trjconv -center would center the COM when it actually
centers the value of (max-min)/2:
http://www.mail-archive.com/gmx-users@gromacs.org/msg41681.html
The short of it is that there are no out-of-the-box tools to construct
the density profile along the normal to a bilayer when the z-dimension
can fluctuate.
You can probably use the trjconv -center (COM) patch after resetting
the box in the following scheme. I used this to create input for
g_spatial, but I suspect that the distribution version of g_density
would also work fine on these files. Note that the only non-standard
thing about this processing is that I applied the trjconv patch to
which I referred above. The idea is to make the box larger than it was
in any of the frames but of a constant size and then center everything
taking careful control of pbc.
rm -rf TEMPORARY_FILES
mkdir -p TEMPORARY_FILES
echo KSC_DOPC | trjconv -f bothsides_center_adjusted_*.xtc -o
/dev/shm/tmp.xtc -n cn.ndx
GMXLIB=/project/pomes/cneale/GPC/exe/intel/gromacs-4.5.3_modtrjconv/exec/share/gromacs/top
echo NE_CZ_NH1_NH2_CB_CG_CD System |
/project/pomes/cneale/GPC/exe/intel/gromacs-4.5.3_modtrjconv/exec/bin/trjconv
-center -pbc mol -f /dev/shm/tmp.xtc -o /dev/shm/tmp2.xtc -s
../../useful/dry.tpr -n cn.ndx
mv /dev/shm/tmp.xtc TEMPORARY_FILES
echo DOPC System |
/project/pomes/cneale/GPC/exe/intel/gromacs-4.5.3_modtrjconv/exec/bin/trjconv
-center -pbc atom -f /dev/shm/tmp2.xtc -o /dev/shm/tmp3.xtc -s
../../useful/dry.tpr -n cn.ndx
mv /dev/shm/tmp2.xtc TEMPORARY_FILES
## now make a new .tpr file in which the solute is at the center of the box
#first output a single frame
echo NE_CZ_NH1_NH2_CB_CG_CD System |
/project/pomes/cneale/GPC/exe/intel/gromacs-4.5.3_modtrjconv/exec/bin/trjconv
-f /dev/shm/tmp3.xtc -dump 125000 -o /dev/shm/tmpgro.gro -s
../../useful/dry.tpr -center -pbc mol -n cn.ndx
#make a new .tpr file
touch empty.mdp
/project/pomes/cneale/GPC/exe/intel/gromacs-4.5.3_modtrjconv/exec/bin/grompp
-p
/project/pomes/cneale/GPC/fromScratch/ARG/MANY_RUNS/TEMPLATE/FILES/complete_dry.top -c /dev/shm/tmpgro.gro -f empty.mdp -o centered.tpr -maxwarn
1
echo NE_CZ_NH1_NH2_CB_CG_CD System |
/project/pomes/cneale/GPC/exe/intel/gromacs-4.5.3_modtrjconv/exec/bin/trjconv
-fit transxy -pbc atom -f /dev/shm/tmp3.xtc -o /dev/shm/tmp4.xtc -s
centered.tpr -n cn.ndx
mv /dev/shm/tmp3.xtc TEMPORARY_FILES
echo System |
/project/pomes/cneale/GPC/exe/intel/gromacs-4.5.3_modtrjconv/exec/bin/trjconv
-box 6 6 14 -f /dev/shm/tmp4.xtc -o /dev/shm/tmp5.xtc -s centered.tpr
-n cn.ndx
mv /dev/shm/tmp4.xtc TEMPORARY_FILES
echo DOPC System |
/project/pomes/cneale/GPC/exe/intel/gromacs-4.5.3_modtrjconv/exec/bin/trjconv
-center -pbc none -f /dev/shm/tmp5.xtc -o /dev/shm/tmp6.xtc -s
centered.tpr -n cn.ndx
mv /dev/shm/tmp5.xtc TEMPORARY_FILES
echo NE_CZ_NH1_NH2_CB_CG_CD System |
/project/pomes/cneale/GPC/exe/intel/gromacs-4.5.3_modtrjconv/exec/bin/trjconv
-fit transxy -pbc none -f /dev/shm/tmp6.xtc -o /dev/shm/tmp7.xtc -s
centered.tpr -n cn.ndx
mv /dev/shm/tmp6.xtc TEMPORARY_FILES
# Now center the solute at 0 0 0
echo NE_CZ_NH1_NH2_CB_CG_CD |
/project/pomes/cneale/GPC/exe/intel/gromacs-4.5.3_modtrjconv/exec/bin/g_traj
-f /dev/shm/tmp7.xtc -ox -s ../../useful/dry.tpr -n cn.ndx -com
avg=($(cat coord.xvg|grep -v '[@|#]' | awk '{x+=$2;y+=$3;z+=$4;n++}
END{print -1*x/n,-1*y/n}'))
/project/pomes/cneale/GPC/exe/intel/gromacs-4.5.3_modtrjconv/exec/bin/trjconv
-f /dev/shm/tmp7.xtc -o /dev/shm/tmp8.xtc -trans ${avg[0]} ${avg[1]} 0
mv /dev/shm/tmp7.xtc TEMPORARY_FILES
Chris.
-- original message --
Hi,
I am trying to calculate the density profile of head group of
bilayer normal to z direction using gromacs 4.0.7. I was trying to
center the density profile about dx/2.dy/2,0 . But, I am finding
problem with using center option. I find using -center option does
not shift the bilayer to 0. The following was my command-lines:
g_density_4mpi -f ../traj_npt -s ../topol -noxvgr -n ../index -b
60 -dens number -center -o density_phosphate_symm.xvg
Any help on how to use the center option will be really helpful.
For this purpose, trjconv is more reliable.
-Justin
--
Justin A. Lemkul
Ph.D. Candidate
ICTAS Doctoral Scholar
MILES-IGERT Trainee
Department of Biochemistry
Virginia Tech
Blacksburg, VA
jalemkul[at]vt.edu | (540) 231-9080
http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin
--
gmx-users mailing listgmx-users@gromacs.org
http://lists.gromacs.org/mailman/listinfo/gmx-users
Please search the archive at
[gmx-users] wibke.sudholt spam needs to stop
That's 9 times in the last 20 days that a single user has hit the list with an auto-reply that seems a lot like an advertisement. Can somebody please ban this user? The most recent example: http://lists.gromacs.org/pipermail/gmx-users/2011-August/063991.html -- 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] more than 100% CPU
Slightly off-topic, but for all but the smallest systems, I get a further 10% efficiency by running a 16-process mpi job on an 8-core machine. I suspect that the story is the same with threads. Thus, on a 16-core node, you could try starting 32 threads. Note that this will report a 2x larger efficiency as you were discussing before, but by checking the ns.day you will see that the benefit is between 5% and 16% Chris. It means you are running a lot of parallel processes, but that does not translate into a linear increase in speed. So faster, but not 16X. Warren Gallin On 2011-08-18, at 5:36 PM, Park, Jae Hyun nmn wrote: Thank you, Warren. Does that mean 16 times faster ? Jae H. Park -Original Message- From: gmx-users-bounces at gromacs.org [mailto:gmx-users-bounces at gromacs.org] On Behalf Of Warren Gallin Sent: Thursday, August 18, 2011 7:26 PM To: Discussion list for GROMACS users Subject: Re: [gmx-users] more than 100% CPU I believe that the current version of GROMACS supports threading, which does not require mpi. So mdrun is running threads at 100% of the activity of each of your 16 nodes, hence 1600%. Warren Gallin On 2011-08-18, at 5:19 PM, Park, Jae Hyun nmn wrote: Hi GMX users, I installed GMX 4.5.3 recently. But, when I just execute mdrun (without mpi, I did not installed mpi-version of mdrun), the use of CPU appears more than 100% (top command in LINUX). How is it possible? For example, I am using 16-node machine. And if I simply run mdrun, then use of CPU is 1600%!!. The simulation runs well and the results looks reasonable. Is there anybody who can teach me what is happening? I would deeply appreciate. -- 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] GROMACS 4.5.4 keep crashing all the time.
You'll need to provide a much better report than this if you want to receive any useful help. Copy and paste the exact commands of what you did Copy and paste the exact log file and error messages Do this for 4.0.7 and 4.5.4, for which I trust that you have been using exactly identical test systems. If not, then please try it again while conserving the system. Chris. -- original message -- Hi Matthew, Thanks for the reply. First, I don't thinks it is poorly compiled bin, as I used both GROMACS compiled by me (on my machine) or by the sys-admin on Linux cluster or blue-gene. The simulations using 4.5.4 crashed giving LINCS error, which is not the case with 4.0.7 (using the same mdp and pdb files). Second, I don't thinks it is poorly compiled bin, as I used both GROMACS compiled with me (on my machine) or by the sys-admin on Linux cluster/blue-gene. cheers, Itamar On 18 August 2011 01:48, Matthew Zwier mczwier at gmail.com wrote: Could be a system blowing up, or perhaps a mis-compiled binary. What error messages do you get when the crash occurs? On Tue, Aug 16, 2011 at 9:48 PM, Itamar Kass itamar.kass at monash.edu wrote: Hi all GROMACS useres and developers, I am interesting in simulating a small protein (~140 aa) in water, with and without Ca ions. In order to do so, I had used version 4.5.4. I had solvate the protein in water, add ions to naturalise the systems, equilibrated the systems and then tried productive runs. Now, no matter what I did, it crashed after few ps's of free MD or during the PR runs. Few of the things I had tried are: 1. Running the simulations on different systems (OSX, linux or blue-gene). 2. Using single or double precision versions. 3. An equilibration stage, 1ns long with a time-step of 1fs, during which the restrained forces where gradually reduced from 1000 to 50 kJ mol-1 nm-2. 4. Running part of the equilibration stage as NVT and then switched to NPT. 5. Started from different x-ray structures, with resolutions differ from 2.5 to 1.7 Angstrom. Finally I had moved back to 4.0.7 which worked like charm. I wonder if someone else had encounter something like this. Attached please find the mdp files I used. All the best, Itamar. -- 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] GROMACS 4.5.4 keep crashing all the time
run an EM with flexible water. I often find that this is the only way to get a stable system. 500 steps of steep with define=-DFLEXIBLE (or different depending on your water model I think) should be enough. Chris. Hi Chris and Justin, On 18/08/2011, at 9:36 AM, Justin A. Lemkul wrote: chris.neale at utoronto.ca wrote: You'll need to provide a much better report than this if you want to receive any useful help. Copy and paste the exact commands of what you did Copy and paste the exact log file and error messages The command I had used (for both 4.0.7 and 4.5.4) are: pdb2gmx -f 1EXR.pdb -o 1EXR.gro -water spc -ignh -p 1EXR.top -i 1EXR.itp editconf -f 1EXR.gro -o 1EXR_box.gro -c -d 1.4 -bt cubic genbox -cp 1EXR_box.gro -cs spc216.gro -o 1EXR_solv.gro -p 1EXR.top -seed 5 grompp -f em.mdp -c 1EXR_solv.gro -p 1EXR.top -o ions.tpr genion -s ions.tpr -o 1EXR_solv_ions.gro -p 1EXR.top -pname NA+ -nname CL- -np 17 -seed 66 grompp -f em.mdp -c 1EXR_solv_ions.gro -p 1EXR.top -o em.tpr mpirun mdrun_d_mpi -v -s em.tpr -x 1EXR_em.xtc -e 1EXR_em.edr -g 1EXR_em.log -c 1EXR_em.gro grompp -f pr.mdp -c 1EXR_em.gro -p 1EXR.top -o pr.tpr mpirun mdrun_d_mpi -v -stepout 1000 -s pr.tpr -x 1EXR_pr.xtc -e 1EXR_pr.edr -g 1EXR_pr.log -o 1EXR_pr.trr -c 1EXR_pr.gro grompp -f md_start.mdp -c 1EXR_pr.gro -p 1EXR.top -o md.tpr mpirun mdrun_d_mpi -v -stepout 1 -s md.tpr -x 1EXR_md.xtc -e 1EXR_md.edr -g 1EXR_md.log -o 1EXR_md.trr -c 1EXR_md.gro The em.mdp: integrator = steep maximum number of steps to integrate nsteps = 5 ;Energy minimizing stuff emstep = 0.001 emtol= 100.0 ; OPTIONS FOR BONDS constraints = none ; OUTPUT CONTROL OPTIONS ; Output frequency for coords (x), velocities (v) and forces (f) nstxout = 1000 nstvout = 1000 nstfout = 1000 ; Output frequency and precision for xtc file nstxtcout= 1000 xtc-precision= 1000 ; Energy monitoring energygrps = system ; NEIGHBORSEARCHING PARAMETERS ; nblist update frequency nstlist = 5 ; ns algorithm (simple or grid) ns-type = Grid ; nblist cut-off rlist= 0.8 ; OPTIONS FOR ELECTROSTATICS AND VDW ; Method for doing electrostatics coulombtype = Reaction-Field rcoulomb = 1.4 ; Dielectric constant (DC) for cut-off or DC of reaction field epsilon-rf = 62 ; Method for doing Van der Waals vdw-type = Cut-off ; cut-off lengths rvdw = 1.4 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = No ; Center of mass control nstcomm = 1000 ; Periodic boundary conditions pbc = xyz ; Mode for center of mass motion removal comm-mode= linear ; Groups for center of mass motion removal comm-grps= system the pr.mdp: define = -DPOSRES integrator = md dt = 0.002 ; ps ! nsteps = 5 ; total 100ps. ; OPTIONS FOR BONDS ; Constrain control constraints = all-bonds ; Do not constrain the start configuration continuation = no ; Type of constraint algorithm constraint-algorithm = lincs ; OUTPUT CONTROL OPTIONS ; Output frequency for coords (x), velocities (v) and forces (f) nstxout = 10 nstvout = 10 nstfout = 10 ; Output frequency and precision for xtc file nstxtcout= 5000 xtc-precision= 1000 ; Energy monitoring energygrps = Protein Non-protein nstenergy= 5000 ; NEIGHBORSEARCHING PARAMETERS ; nblist update frequency nstlist = 5 ; ns algorithm (simple or grid) ns-type = Grid ; nblist cut-off rlist= 0.8 ; OPTIONS FOR ELECTROSTATICS AND VDW ; Method for doing electrostatics coulombtype = Reaction-Field rcoulomb = 1.4 epsilon_rf = 62; As suggested at J Comput Chem. 2004 Oct;25(13):1656-76. vdw-type = Cut-off ; cut-off lengths rvdw = 1.4 ; Apply long range dispersion corrections for Energy and Pressure DispCorr = no ; OPTIONS FOR WEAK COUPLING ALGORITHMS ; Temperature coupling tcoupl = berendsen ; Groups to couple separately, time constant (ps) and reference temperature (K) tc-grps = Protein Non-Protein tau-t= 0.1 0.1 ref-t= 300 300 ; Pressure coupling Pcoupl = berendsen Pcoupltype = isotropic ; Time constant (ps), compressibility (1/bar) and reference P (bar) tau_p= 1.0 compressibility = 4.5e-5 ref_p= 1.0 ; Generate velocites is on at 290K
[gmx-users] GROMACS 4.5.4 keep crashing all the time
OK, here's my last few ideas: 1. Please try to repeat this with gen_vel set to the same value as your temperature coupling 2. Can you reproduce this in serial? 3. Can you reproduce this with the sd integrator? 4. Can you reproduce this with a simpler system? protein in vacuum or just water or remove the ions, etc? 5. Take the output .gro from 4.0.7 that ran fine for X ns and run it under 4.5.4. Do you get the same lincs warnings? 6. Also, note that you are getting warnings and the run does not actually crash but just stops after too many warnings. So what are atoms 981 and 982? Does their motion look different in an important ways between the 4.0.7 and 4.5.4 trajectories? Chris. -- original message -- Hi Chris, thanks for the advice, I have to say I tried this as well without any success. Itamar On 18/08/2011, at 11:11 AM, chris.neale at utoronto.ca wrote: run an EM with flexible water. I often find that this is the only way to get a stable system. 500 steps of steep with define=-DFLEXIBLE (or different depending on your water model I think) should be enough. Chris. Hi Chris and Justin, On 18/08/2011, at 9:36 AM, Justin A. Lemkul wrote: -- 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] suggestion that mdrun should ensure npme the number of processes
Currently, gromacs4.5.4 gives a segfault if one runs mpirun -np 8 mdrun_mpi -npme 120 with no warning of the source of the problem. Obviously npmennodes is a bad setup, but a check would be nice. Thank you, Chris. -- 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] Constraints not working in pull code (sometimes, sometimes not)
Dear Krapnik: 1. please make the test cases identical, that means doing a simulation that you may not really be interested in so that the solutes are the same for both lipid and octane. It also means setting the displacements to similar values in my opinion (because perhaps the problem is based on your box size along z and choice of pull_pbcatom for example). There is *something* strange happening, but let's eliminate all other possibilities before assuming that there is a problem with the source code. 2. Once you have a setup in which the only difference is octane vs. lipid (really the only difference please -- even add your octane pressure settings to the lipid simulation) and you find strange behaviour for octane and not lipid, then please attempt that same system with pull=umbrella in place of pull=constraints and see if the solute also drifts far from its initial position. 3. Your initial definition of the problem was a little misleading. Both Justin and I focused on third decimal place differences, but now you are saying that there is nm-scale motion along z. Please report those values for your test systems, which are outlined above, after you redo your tests with identical settings other than the switch from lipid to octane. -- I no longer thing that a double precision test is necessary if you are seeing nm-scale movement along z. 4. Perhaps I was a little hasty to complain about a bug-report style post. I did like all of the information that you provided. I just got the impression from your title and the end of your post that you thought it was a problem with gromacs (which it may still be) and I think that this is generally a counter-productive attitude because it hinders one from making the proper tests. Chris -- original message -- Re: Constraints not working in pull code (sometimes, sometimes not) Krapnik krapnik at gmail.com Thu Aug 11 09:53:24 CEST 2011 * Previous message: [gmx-users] append files in Gromacs 3.3 * Next message: [gmx-users] Constraints not working in pull code (sometimes, sometimes not) * Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] Dear Justin Thank you for your responses, but it still is not clear to me why is that happening I don't see any real problem here. The DOPC membrane probably holds its shape better than the octanol slab and thus the reference position (which determines the constraint) is less mobile. That is certain. When we were running free simulation, the box started to narrow. Therefore we have frozen the box in xy plane. As such, the constraint is more stable in DOPC than in octanol. The slight bump in the octanol system amounts to a change of 0.16%, which I'd say is quite small. The trouble is, that during next 10 ns the distance was between 2.0 - 4.0 and molecule moved outside of the ocatnol slab and inside, which is then useless for calculation of PMF. You're also looking at the first few frames of the simulation, which may amount to equilibration, but you haven't mentioned if you've done anything prior to this run. Previously, a free simulation with molecule on slab was run for 20 ns, from which frame with small molecule at specidied depth was taken, so I suppose, that the start is equilibrated. As long as there is no systematic change in position, I'd say there's no problem. Unfortunately, there is as I have said before. -Justin Dear Chris, thank you for responses as well. I agree that Justin is probably correct, although constraints should technically work just fine with a highly dynamic reference group. Any problems should show up as the system blowing up, but not in correctly setting the position. I hope, that the position is set up ok, as it stays steady in DOPC. I think that the problems can arise, however, when you have multiple competing constraints. You might, for example, try without constraining the octane bonds -- I think that you should get perfect match in that case. I will try. The erroneous movement of molecule seems to be enhanced within XY-plane frozen. The curious thing is however, that when I did not constrain size of xy-plane, then the distance to the COM of slab is not changing much (only last digit changes), which in comparison to frozen XY-plane, where there is huge movement above the slab and within is rather ok to me. Unfortunately, I have the issue with narrowing of the box in the unconstrained XYplane simulation, which also destroys any attempts for PMF as I have mentioned before. Also, it is worth comparing in single and double precision. I will try this either. Still, I am not sure that you actually did a proper comparison for the following 2 reasons: 1. you have apparently different masses for the pulled group: Pull group 1:15 atoms, mass 194.194 Pull group 1:15 atoms, mass 146.146 The difference is not much in my opinion, while molecules tested are similar in size and are run
[gmx-users] vsites with NAC and opls
This is fixed, just posting for posterity. If one wants to run pdb2gmx with -vsites hydrogens on gromacs 4.5.4 or 4.0.7 while using the oplss ff, there is the error message: Fatal error: Can't find dummy mass for type opls_242 bonded to type opls_238 in the virtual site database (.vsd files). Add it to the database! To fix this, one must then add the following line: opls_242 opls_238 MCH3A to the [ CH3 ] section of oplsaa.ff/aminoacids.vsd -- 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] Constraints not working in pull code (sometimes, sometimes not)
I agree that Justin is probably correct, although constraints should technically work just fine with a highly dynamic reference group. Any problems should show up as the system blowing up, but not in correctly setting the position. I think that the problems can arise, however, when you have multiple competing constraints. You might, for example, try without constraining the octane bonds -- I think that you should get perfect match in that case. Also, it is worth comparing in single and double precision. Still, I am not sure that you actually did a proper comparison for the following 2 reasons: 1. you have apparently different masses for the pulled group: Pull group 1:15 atoms, mass 194.194 Pull group 1:15 atoms, mass 146.146 2. you used very different starting depth offsets: -2.81855 is not very close to -1.60929 PS: I would personally prefer if you avoided jumping directly to a bug report unless you are sure of it. There is always at least the possibility that you are not using the code correctly. Chris. Dear all, we are trying to simulate molecule in specified depth in octanol to calculate logP. for this purpose we have used box with slab of octanol, molecule in specified distance in z-axis and this mdp: integrator = md dt = 0.002 tinit = 0 nsteps = 5125000 ; 10.250 ns nstcomm = 1 comm_mode = Linear ; Output parameters nstxout = 0 ; every 1 ns nstvout = 0 nstfout = 0 nstxtcout = 500 ; every 1 ps nstenergy = 500 ; Bond parameters constraints = all-bonds ; Single-range cutoff scheme nstlist = 5 ns_type = grid rlist = 1.4 rcoulomb= 1.4 rvdw= 1.4 ; PME electrostatics parameters coulombtype = PME fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-5 optimize_fft= yes ; Berendsen temperature coupling is on in two groups Tcoupl = V-rescale tc_grps = System tau_t = 0.1 ref_t = 310 ; Pressure coupling is on Pcoupl = berendsen pcoupltype = anisotropic ; anisotropic pressure coupling tau_p = 10 compressibility = 0 0 4.5e-5 0 0 0 ; since octanol tries to separate from water and shortening xy plane ref_p = 1.0 1.0 1.0 0 0 0 ; Pull code pull= constraint; pull_geometry = distance ; Pull along the vector connecting the two groups. Components can be selected with pull_dim. pull_dim= N N Y ; put potential only to axis z pull_start = yes ; define initial COM distance 0 pull_ngroups= 1 ; one group to pull pull_group0 = OCT ; reference group (can be empty to set it to 0,0,0) pull_group1 = DRG pull_rate1 = 0 ; 0.01 nm per ps = 10 nm per ns = 10 m per s pull_k1 = 2000 ; kJ mol^-1 nm^-2 however when I list pull-x.xvg I obtain: @ s0 legend 0 Z @ s1 legend 1 dZ 0. 6.06871 -2.81855 0.0200 6.07793 -2.81861 0.0400 6.08787 -2.81856 0.0600 6.08462 -2.81855 0.0800 6.10028 -2.82007 The situation is even more strange when almost the same mdp with DOPC membrane works fine as can be seen in pull-x.xvg: 0. 3.72735 -1.60929 0.0200 3.72727 -1.60929 0.0400 3.7272 -1.60929 0.0600 3.72716 -1.60929 0.0800 3.72715 -1.60929 0.1000 3.72717 -1.60929 Differences in the working case in mdp are only this: compressibility = 4.5e-5 4.5e-5 4.5e-5 0 0 0 ref_p = 1.0 1.0 1.0 0 0 0 pull_group0 = DOPC Strangerly enough - the same happens in GMX4.0.7 and in GMX4.5.3 I wonder where the difference is, as both logs states that Will apply constraint COM pulling in geometry 'distance' between a reference group and 1 group Pull group 0: 6912 atoms, mass 100624.859 Pull group 1:15 atoms, mass 194.194 Will apply constraint COM pulling in geometry 'distance' between a reference group and 1 group Pull group 0: 9570 atoms, mass 124631.453 Pull group 1:15 atoms, mass 146.146 Where should I look and repair? Thank you in advance Karel -- 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] spatial distribution function in a binary solvent mixture
I think that you have a misconception about what g_spatial does. For a system with many type A and many type B, you need to average over all of one type as the central solute to compute an rdf, and perhaps that is what you want for your sdf. g_spatial, however, does not do any fitting. g_sdf did that. You can obtain an old version of the source (4.0.7 should have it I think) if you want to use g_sdf. I have never used g_sdf myself. In case that doesn't answer your question, then let me make one more point: g_spatial requires that a bin exists for every count. Thus either (a) find a large memory node somewhere or (b) pre-process your trajectory using trjorder to order the solvent molecules and then keep only the N closest, then run g_spatial. But again, I think that this will not provide what you are looking for but is likely to give you a smear over your box. Chris. -- original message -- Dear all, I'm working with a binary solvent mixture containing 2000 molecules (1800 type-A + 200 type-B). Both the types of solvent molecules have similar structure (they are both diatomic molecules) except the polarity.I'm trying to calculate sdf of type-B solvent molecules. I followed the step-by-step instructions from the manual using g_spatial. I'm trying to reduce the bin width to 0.05A. The *.cube file generated with a bin width of 0.09A is already 4.9GB in size. As I'm more interested in the first solvation shell around the type-B solvent molecule, I was wondering if I could find a way to control the maximum radius of the sphere around central molecule (center of coordinates?)? Also, can anyone please let me know how g_spatial deals with the angular part of sdf? (I've searched a lot but could not gather enough meaningful information) regards, Debasmita -- 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] spatial distribution function in a binary solvent mixture
Here's the totally wasteful way for you to get what you want: 1. for each molecule in type A: trjorder everything around that molecule and output only the closest N atoms of type A and the closest M atoms of type B. Ensure that the coordinate order is: your central molecule of type A is first, the remaining type A second, and the type B third. 2. concatenate all of these trajectories into one single large trajectory 3. trjconv -fit trans+rot , fitting only to the single central molecule. 4a. g_spatial for the non-central type A 4b. g_spatial for type B * now repeat switching types A and B in the above steps to get the SDFs arouns type B. A much better idea is for you to modify g_spatial to do this all within one loop. But if you don't know C then you are out of luck there. g_sdf may not work because you are underdetermined (there is a symmetry about the long axis that you can not specify). But perhaps you could use that to your advantage and select some totally unrelated atom as the third atom for the fit. My best guess is that this would be no worse than g_spatial. A final option is for you to use g_mindist to output all of the contacts and then run over that with your own script to process the data into an SDF and output it in cube format. Chris. -- original message -- Thanks Chris. I presume g_sdf won't be helpful for my system. [Hide Quoted Text] I think that you have a misconception about what g_spatial does. For a system with many type A and many type B, you need to average over all of one type as the central solute to compute an rdf, and perhaps that is what you want for your sdf. g_spatial, however, does not do any fitting. g_sdf did that. You can obtain an old version of the source (4.0.7 should have it I think) if you want to use g_sdf. I have never used g_sdf myself. In case that doesn't answer your question, then let me make one more point: g_spatial requires that a bin exists for every count. Thus either (a) find a large memory node somewhere or (b) pre-process your trajectory using trjorder to order the solvent molecules and then keep only the N closest, then run g_spatial. But again, I think that this will not provide what you are looking for but is likely to give you a smear over your box. Chris. -- original message -- Dear all, I'm working with a binary solvent mixture containing 2000 molecules (1800 type-A + 200 type-B). Both the types of solvent molecules have similar structure (they are both diatomic molecules) except the polarity.I'm trying to calculate sdf of type-B solvent molecules. I followed the step-by-step instructions from the manual using g_spatial. I'm trying to reduce the bin width to 0.05A. The *.cube file generated with a bin width of 0.09A is already 4.9GB in size. As I'm more interested in the first solvation shell around the type-B solvent molecule, I was wondering if I could find a way to control the maximum radius of the sphere around central molecule (center of coordinates?)? Also, can anyone please let me know how g_spatial deals with the angular part of sdf? (I've searched a lot but could not gather enough meaningful information) regards, Debasmita -- 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] umbrella sampling
I've never actually reweighted with F-values. These are the values by which the component PMFs are shifted during WHAM after debiasing the umbrella potentials. Look at equations 6 and 7 in The multicanonical weighted histogram analysis method for the free-energy landscape along structural transition paths Satoshi Ono, Nobuyuki Nakajima, Junichi Higo and Haruki Nakamura My best guess is that you do it like this: 1. Run 1-D wham and obtain the F-values. 2. create an empty 3D histogram structure for z, rgyrA, and rgyrB. 3. Go through the data from each window and add a value to the appropriate location for [z,rgyrA,rgyrB]. This value is not the integer count 1. This value that you add to each histogram bin is a real number: 1 * exp(-(1/RT)*-F) * exp(-(1/RT)*0.5K*(z-z0)^2 4. Now convert your histogram of probabilities to free energies dG=-RTln(P) ** To test this, project your 3-D free energy profile onto each dimension successively and compare it to (z) the profile from 1-D wham, and (rgyrA or rgyrB) the profile from 2-D wham with zero force constants as I outlined in my previous post. PLEASE NOTE: this email outlines an *idea* of how to accomplish what you want. You must check the match for yourself. Chris. Qian Wang qwang at mail.uh.edu Sat Jul 30 18:05:33 CEST 2011 * Previous message: [gmx-users] umbrella sampling * Next message: [gmx-users] OPLSAA parameters * Messages sorted by: [ date ] [ thread ] [ subject ] [ author ] Thanks very much. I am sorry I did not describe my question very clear. In my system, I am going to obtain three things. (1) the free energy versus distance, which I can obtain from g_wham. (2) the free energy versus, eg, the radius of gyration of the object A. (3) the 2-D free energy as a function of the radius of gyration of A and the radius of gyration of B. I do not know how to get (2) and (3). Your first method can help me obtain (2), but I can not obtain (3) because then it needs a 3-D WHAM code, am I correct? For your second method, I do not quite understand how to re-weighting using the F-values from 1D-WHAM. Could you please explain it? Thanks. Sincerely, Qian - Original Message - From: chris.neale at utoronto.ca Date: Friday, July 29, 2011 11:00 pm Subject: [gmx-users] umbrella sampling To: gmx-users at gromacs.org format your data for 2D-WHAM with 1D being the distance and the 2nd-D being your other coordinate of interest. Specify a value of zero for the force constants for your 2nd-D. Run 2D-WHAM. Boltzmann project the 2D PMF onto your 2nd-D. I think you can also do essentially the same thing by re- weighting using the F-values from 1D-WHAM, but I find the above method to be the simplest. It also provides you with a 2D free energy profile, which can be informative both biologically and to indicate on sampling problems. Note that you're very likely going to run into convergence problems since your 2nd-D will rely on brute-force to converge, and worse: the umbrellas in 1D can force the sampling in the 2nd- D to surmount energy barriers that might be circumvented in unrestrained sampling. Chris. -- original message -- Qian Wang wrote: -- 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] umbrella sampling
format your data for 2D-WHAM with 1D being the distance and the 2nd-D being your other coordinate of interest. Specify a value of zero for the force constants for your 2nd-D. Run 2D-WHAM. Boltzmann project the 2D PMF onto your 2nd-D. I think you can also do essentially the same thing by re-weighting using the F-values from 1D-WHAM, but I find the above method to be the simplest. It also provides you with a 2D free energy profile, which can be informative both biologically and to indicate on sampling problems. Note that you're very likely going to run into convergence problems since your 2nd-D will rely on brute-force to converge, and worse: the umbrellas in 1D can force the sampling in the 2nd-D to surmount energy barriers that might be circumvented in unrestrained sampling. Chris. -- original message -- Qian Wang wrote: Hi, I used umbrella sampling method to restrain the distance of two molecules at several distances. Then I can use g_wham to get the free energy as a function of the distance. Is there any way that I can get the free energy as a function of another parameter? Thanks a lot. -- 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] large error bars in PMF
I don't see why the box-type makes any difference whatsoever. It is possible that if you use a rhombic dodecahedron, you may reduce the system size, thus simulate more ns/day, thus converge faster, but that should be the only effect. I would be interested to hear more from Justin about how the box-shape is expected to affect peptide rotation... perhaps I misunderstand this point. I have a few other comments. 1. If you allow the peptide to rotate freely, then you do indeed need to converge all of their different rotational interactions. An alternative is to apply orientational restraints during the pulling (assuming that you know the bound state) and then to correct to an unrestrained state at large separations. You can see, for instance, D. L. Mobley, J. D. Chodera, K. A. Dill. On the use of orientational restraints and symmetry number corrections in alchemical free energy calculations, ... 2. You are using pull_dim = N N Y which, if I haven't entirely forgotten how the pull-code works, means that the distance along Z is restrained but the distance along X and Y is free to change. What you end up with by sampling in this way is pretty strange and will require a really weird volumetric correction in the absence of infinite sampling time. You must decide to either: (i) pull to a spherical distance: pull_dim = Y Y Y pull_geometry = distance or (ii) to pull along a defined vector pull_dim = Y Y Y pull_geometry = direction What you have done: pull_dim = N N Y pull_geometry = distance is only really useful when the system is isotropic along the XY plane (at least in the time averaged sense), such as for a lipid bilayer, or when the freedom in X and Y is very low, such as in a channel. 3. Finally, just because you sampled *more* doesn't mean that your values are converged. Look into block averaging and test to see if your binding free energy is drifting over time. Good luck, Chris. -- original message -- Hi again, I have one doub about the suggestion of using a dodecahedral box for my umbrella sampling to remove the problems I am having with the peptides rotating. I cannot see why a dodecaheral box is going to avoid this. Would it be better a truncated octahedron? Thanks a lot for your help. Best wishes, Rebeca. [Hide Quoted Text] Date: Thu, 21 Jul 2011 15:16:52 -0400 From: jalem...@vt.edu To: gmx-users@gromacs.org Subject: Re: FW: [gmx-users] large error bars in PMF Rebeca García Fandiño wrote: I am trying to achieve the binding energy of the dimer composed by the two small cyclic peptides, to compare it with experimental. What advantages would I have using 3D PMF instead only 1D for this calculation? Intuitively, two molecules diffuse through solution until they find one another, which to me sounds a lot like a 3D path. Further, using a dodecahedral box for your umbrella sampling removes the problems you're having with the peptides rotating. It sounds like you're trying to pull in one direction along a rectangular box, but the peptides are not playing nice. I feel like this discussion has come up at least once or twice before, though... -Justin Thanks a lot! Rebeca. Date: Thu, 21 Jul 2011 14:14:44 -0400 From: jalem...@vt.edu To: gmx-users@gromacs.org Subject: Re: [gmx-users] large error bars in PMF Rebeca García Fandiño wrote: Hi, thanks a lot for your quick answer. What I am trying to pull are two small peptides one from another (r_1 and r_2). I did not understand very well your last suggestion: ...if you want reasonable error bars you will not lots of well-converged data. Oops, that should have read you will *need* lots of well-converged data. Do you mean I will need also more windows besides extending the simulations? I doubt you need more windows. Likely you just need more time in each. I think the problem could be also that the peptides I am using rotate in the box and they do not remain flat one respect to the other. They gyrate freely and some parts of their structure interact along the pulling... Interactions are part of the dissociation process and are not problematic per se. But if you're trying to obtain only a one-dimensional PMF then your rotation could be a problem. Is there some reason you need a one-dimensional PMF and not a three-dimensional PMF? What are you trying to achieve? -Justin Thanks a lot again for your help. Best wishes, Rebeca. From: rega...@hotmail.com To: gmx-users@gromacs.org Date: Thu, 21 Jul 2011 16:36:59 + Subject: [gmx-users] large error bars in PMF Hi, I am trying to calculate the binding energy of two molecules using the PMF (Umbrella Sampling method) and Gromacs 4.0. Some weeks ago I have written to the list because changing the number
[gmx-users] large error bars in PMF
I see now what you mean. As it happens, I doubt that this would have caused the problem since no force was applied on X and Y dimensions, so it would require that there was a PBC-based distance degeneracy along Z, although this is of course possible and hopefully Rebecca will answer this part. Also, thanks for the pull_dimension/pull_vec fix. Chris. -- original message -- chris.ne...@utoronto.ca wrote: [Hide Quoted Text] I don't see why the box-type makes any difference whatsoever. It is possible that if you use a rhombic dodecahedron, you may reduce the system size, thus simulate more ns/day, thus converge faster, but that should be the only effect. I would be interested to hear more from Justin about how the box-shape is expected to affect peptide rotation... perhaps I misunderstand this point. My point was not that the box shape has any effect on protein rotation. That will happen regardless of the box shape, of course. My suggestion for this box type was that since Rebeca has a system that is essentially spherically symmetric (i.e. two proteins connected by some arbitrary vector, which are at the same time rotating freely), then she must use a suitable box shape that reflects this type of symmetry. I never got a clear answer to whether or not the weird interactions she cited were due to PBC or not, but if one uses a rectangular box for a system like this one, there can be artificial interactions very easily. [Hide Quoted Text] I have a few other comments. 1. If you allow the peptide to rotate freely, then you do indeed need to converge all of their different rotational interactions. An alternative is to apply orientational restraints during the pulling (assuming that you know the bound state) and then to correct to an unrestrained state at large separations. You can see, for instance, D. L. Mobley, J. D. Chodera, K. A. Dill. On the use of orientational restraints and symmetry number corrections in alchemical free energy calculations, ... 2. You are using pull_dim = N N Y which, if I haven't entirely forgotten how the pull-code works, means that the distance along Z is restrained but the distance along X and Y is free to change. What you end up with by sampling in this way is pretty strange and will require a really weird volumetric correction in the absence of infinite sampling time. You must decide to either: (i) pull to a spherical distance: pull_dim = Y Y Y pull_geometry = distance or (ii) to pull along a defined vector pull_dim = Y Y Y pull_geometry = direction Just a note here - if you set direction geometry, the pull_dim keyword is not used, but pull_vec is. -Justin [Hide Quoted Text] What you have done: pull_dim = N N Y pull_geometry = distance is only really useful when the system is isotropic along the XY plane (at least in the time averaged sense), such as for a lipid bilayer, or when the freedom in X and Y is very low, such as in a channel. 3. Finally, just because you sampled *more* doesn't mean that your values are converged. Look into block averaging and test to see if your binding free energy is drifting over time. Good luck, Chris. -- original message -- Hi again, I have one doub about the suggestion of using a dodecahedral box for my umbrella sampling to remove the problems I am having with the peptides rotating. I cannot see why a dodecaheral box is going to avoid this. Would it be better a truncated octahedron? Thanks a lot for your help. Best wishes, Rebeca. ... snip... -- 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] g_msd bug
600 GB of memory? I highly doubt that you have that much memory available. Are you sure that this is not a typo? Can you please post evidence that you have =600 GB of memory available? It is common for clusters to disallow an individual process from using 10% of the total memory on a head-node, which makes 600 MB more likely, in which case you can try submitting your job to a compute node. -- original message -- Hello, thank you for your reply. I have used following command : g_msd -n POPC.ndx -lateral z -o POPC_msd.xvg -mol POPC_diff.xvg Trajectory has 1 frames and the system it was ran on is Fedora Red Hat 5.4. Indeed my network administrator was very unhappy about comsumed memory. Regards, Slawomir Wiadomo¶æ napisana przez Tsjerk Wassenaar w dniu 2011-06-27, o godz. 15:40: [Hide Quoted Text] Hi Slawomir, That's quite a usage of memory! Can you provide more information? Like the number of frames in the trajectory, the command line you used, and the system you ran on? Cheers, Tsjerk 2011/6/27 S³awomir Stachura stachura.slawo...@gmail.com: Hi GMX Users, I am writting this email, beacause I think the g_msd program in Gromacs 4.5.4 bears a problem. I was calculating the MSD od center of mass of POPC in membrane (system contains 274 POPC lipid molecules in all-atom force field) from 50 ns trajectory and it seems to consume great amount of memory. With time of calculations the memory reserves are gradually devoured to the extent, in my case, of over 600 GB (than my administrator of cluster killed the process). It seems that it does not release memory and it's pilling results up with steps in memory. Have you heard of such case? Best wishes, -- 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] gromacs 2.1
Thank you Rossen and Roland.
I was able to obtain the source but not compile it. It would be
appreciated if you can offer any suggestions.
I obtained gromacs via:
git clone git://git.gromacs.org/gromacs.git
git checkout release-2-0-01
But I can not compile it. When I cd src and run make or gmake I get
the errors listed below.
... snip ...
gmake[348]: Entering directory
`/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src'
cd include; gmake oclean; cd ..; cd gmxlib; gmake oclean; cd ..; cd
mdlib; gmake oclean; cd ..; cd fftw; gmake oclean; cd ..; cd kernel;
gmake oclean; cd ..; cd tools; gmake oclean; cd ..;
/bin/sh: line 0: cd: include: No such file or directory
gmake[349]: Entering directory
`/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src'
cd include; gmake oclean; cd ..; cd gmxlib; gmake oclean; cd ..; cd
mdlib; gmake oclean; cd ..; cd fftw; gmake oclean; cd ..; cd kernel;
gmake oclean; cd ..; cd tools; gmake oclean; cd ..;
... snip ...
I tried the documentation, but the FAQ file at
file:///home/cneale/work/June2011/gromacs/share/html/gmxfaq.html#install links
to a non-existant file
file:///home/cneale/work/June2011/gromacs/share/html/readme2.0.html#install
.
I moved the install directory into src/ and that solved the error
messages that I was obtaining, but now make errors with:
/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src/include/smalloc.h:36:2: warning: #ident is a deprecated GCC
extension
dlist.c:197:1: error: pasting . and H does not give a valid
preprocessing token
dlist.c:197:1: error: pasting . and N does not give a valid
preprocessing token
dlist.c:197:1: error: pasting . and Cn does not give a valid
preprocessing token
dlist.c:197:1: error: pasting . and C does not give a valid
preprocessing token
dlist.c:200:1: error: pasting . and N does not give a valid
preprocessing token
dlist.c:200:1: error: pasting . and Cn does not give a valid
preprocessing token
dlist.c:200:1: error: pasting . and C does not give a valid
preprocessing token
dlist.c:200:1: error: pasting . and O does not give a valid
preprocessing token
dlist.c:203:1: error: pasting . and minO does not give a valid
preprocessing token
dlist.c:203:1: error: pasting . and minC does not give a valid
preprocessing token
dlist.c:203:1: error: pasting . and N does not give a valid
preprocessing token
dlist.c:203:1: error: pasting . and Cn does not give a valid
preprocessing token
dlist.c:212:1: error: pasting . and Cn does not give a valid
preprocessing token
dlist.c:212:1: error: pasting . and Cn does not give a valid
preprocessing token
dlist.c:212:1: error: pasting . and Cn does not give a valid
preprocessing token
dlist.c:212:1: error: pasting . and Cn does not give a valid
preprocessing token
make[1]: *** [dlist.o] Error 1
make[1]: Leaving directory
`/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src/tools'
Also moving the lib directory into the src directory did not help any further.
The compiler was gcc version 4.4.0
Thank you very much,
Chris.
-- original message --
You can try also
$ git tag
release-1-6
release-1-6-01
release-2-0
release-2-0-01
release-2-0-02
release-3-0
release-3-1
release-3-1-1
release-3-1-2
release-3-1-3
release-3-1-4
release-3-2
v4.0.7
v4.5
v4.5-beta1
v4.5-beta2
v4.5-beta3
v4.5-beta4
v4.5.1
v4.5.2
v4.5.3
v4.5.4
# git checkout release-2-0-01
On 6/12/11 11:49 AM, Roland Schulz wrote:
Hi,
you can get it from git. git log {file} shows the history and git
show {version}:{file} gives you a specific version.
In this case you can use:
git show 74b20ce9:src/tools/g_order.c
Is that the official 2.1? We should tag it then.
Rossen
Roland
On Sat, Jun 11, 2011 at 10:11 PM, chris.neale at utoronto.ca
mailto:chris.neale at utoronto.ca wrote:
Dear users:
Does anybody have the source code for gromacs version 2.1? I would
like to check the original source of g_order, but only versions
3.0 and above are available on the gromacs website.
Thank you,
Chris.
-- gmx-users mailing list gmx-users at gromacs.org
mailto:gmx-users at 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 thewww
interface or send it to gmx-users-request at gromacs.org
mailto:gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
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] gromacs 2.1
In the middle of my last post, it should have stated:
I moved the **include** directory into src/ and that solved the error
messages that I was obtaining, but now make errors with:
Sorry for the confusion,
Chris.
Quoting chris.ne...@utoronto.ca:
Thank you Rossen and Roland.
I was able to obtain the source but not compile it. It would be
appreciated if you can offer any suggestions.
I obtained gromacs via:
git clone git://git.gromacs.org/gromacs.git
git checkout release-2-0-01
But I can not compile it. When I cd src and run make or gmake I get
the errors listed below.
... snip ...
gmake[348]: Entering directory
`/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src'
cd include; gmake oclean; cd ..; cd gmxlib; gmake oclean; cd ..; cd
mdlib; gmake oclean; cd ..; cd fftw; gmake oclean; cd ..; cd kernel;
gmake oclean; cd ..; cd tools; gmake oclean; cd ..;
/bin/sh: line 0: cd: include: No such file or directory
gmake[349]: Entering directory
`/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src'
cd include; gmake oclean; cd ..; cd gmxlib; gmake oclean; cd ..; cd
mdlib; gmake oclean; cd ..; cd fftw; gmake oclean; cd ..; cd kernel;
gmake oclean; cd ..; cd tools; gmake oclean; cd ..;
... snip ...
I tried the documentation, but the FAQ file at
file:///home/cneale/work/June2011/gromacs/share/html/gmxfaq.html#install
links
to a non-existant file
file:///home/cneale/work/June2011/gromacs/share/html/readme2.0.html#install
.
I moved the install directory into src/ and that solved the error
messages that I was obtaining, but now make errors with:
/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src/include/smalloc.h:36:2: warning: #ident is a deprecated
GCC
extension
dlist.c:197:1: error: pasting . and H does not give a valid
preprocessing token
dlist.c:197:1: error: pasting . and N does not give a valid
preprocessing token
dlist.c:197:1: error: pasting . and Cn does not give a valid
preprocessing token
dlist.c:197:1: error: pasting . and C does not give a valid
preprocessing token
dlist.c:200:1: error: pasting . and N does not give a valid
preprocessing token
dlist.c:200:1: error: pasting . and Cn does not give a valid
preprocessing token
dlist.c:200:1: error: pasting . and C does not give a valid
preprocessing token
dlist.c:200:1: error: pasting . and O does not give a valid
preprocessing token
dlist.c:203:1: error: pasting . and minO does not give a valid
preprocessing token
dlist.c:203:1: error: pasting . and minC does not give a valid
preprocessing token
dlist.c:203:1: error: pasting . and N does not give a valid
preprocessing token
dlist.c:203:1: error: pasting . and Cn does not give a valid
preprocessing token
dlist.c:212:1: error: pasting . and Cn does not give a valid
preprocessing token
dlist.c:212:1: error: pasting . and Cn does not give a valid
preprocessing token
dlist.c:212:1: error: pasting . and Cn does not give a valid
preprocessing token
dlist.c:212:1: error: pasting . and Cn does not give a valid
preprocessing token
make[1]: *** [dlist.o] Error 1
make[1]: Leaving directory
`/project/pomes/cneale/GPC/exe/intel/gromacs-git/gromacs/src/tools'
Also moving the lib directory into the src directory did not help
any further.
The compiler was gcc version 4.4.0
Thank you very much,
Chris.
-- original message --
You can try also
$ git tag
release-1-6
release-1-6-01
release-2-0
release-2-0-01
release-2-0-02
release-3-0
release-3-1
release-3-1-1
release-3-1-2
release-3-1-3
release-3-1-4
release-3-2
v4.0.7
v4.5
v4.5-beta1
v4.5-beta2
v4.5-beta3
v4.5-beta4
v4.5.1
v4.5.2
v4.5.3
v4.5.4
# git checkout release-2-0-01
On 6/12/11 11:49 AM, Roland Schulz wrote:
Hi,
you can get it from git. git log {file} shows the history and git
show {version}:{file} gives you a specific version.
In this case you can use:
git show 74b20ce9:src/tools/g_order.c
Is that the official 2.1? We should tag it then.
Rossen
Roland
On Sat, Jun 11, 2011 at 10:11 PM, chris.neale at utoronto.ca
mailto:chris.neale at utoronto.ca wrote:
Dear users:
Does anybody have the source code for gromacs version 2.1? I would
like to check the original source of g_order, but only versions
3.0 and above are available on the gromacs website.
Thank you,
Chris.
-- gmx-users mailing list gmx-users at gromacs.org
mailto:gmx-users at 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 thewww
interface or send it to gmx-users-request at gromacs.org
mailto:gmx-users-request at gromacs.org.
Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
--
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
[gmx-users] pulling code
Dear Adam: I have a number of questions listed below. But first, I think it would be useful if you explain exactly what you want to accomplish and how you tried to do this with non-modified code and what problem you ran into that caused you to modify the code. After you outline that, please address the following questions: 1. Why not use one terminus as the reference group and the other terminus as the pull_group0 ? 2. How did you modify the code to use the shortest vector? I was under the impression that the code did already do this (shortest vector accounting for PBC) and that this was actually something you wanted to avoid. 3. Why can you not just use A-B even if A is positive and B is negative? If you are trying to avoid determining a PBC-wraparound distance, then it seems like A-B does give you exactly what you want as long as the entire protein is inside the unit cell. 4. If you really need to know the box size, t is available in the code at runtime. Why can you not just use the appropriate variable? 5. How are you fixing the non-terminal groups that you say need to be fixed? I don't understand why you can't fix groups with position restraints and the pull with the pull code. Also, why do you need to fix certain non-terminal groups? 6. This confused me a lot: I don't think I want to freeze or restrain the groups that need to be fixed, because they must still be allowed flexibility to twist and deform--I only want to roughly fix their centers of mass. ... so then how did you fix them? Chris. -- original message -- Hi all, I am trying to measure the energy curve of unwinding the SNARE protein complex. This requires pulling the C-termini of two helices apart from each other, while fixing certain other parts of the complex. To do this I excluded the reference group, so that all groups were pulled to absolute coordinates. However, as the protein was centered in the box, and because the box center is where the coordinates transition from -boxlength/2 to +boxlength/2 when u7. Why do you need to fix certain non-sing pbc (ie. the dateline), every time a group passes through the center of the box in any of the 3 dimensions, it is suddenly yanked violently through the box because it thinks it is now a full boxlength away from its destination. I modified the pull code to take the shortest vector between the current position and destination, which seems to fix this problem. I also modified g_wham so that it can take the difference between two groups' positions as the energy curve coordinate (since several groups need to be fixed, neither of these two groups can be used as the reference group). However, as I use pressure coupling, the box dimensions change during the simulation, and if one group has a positive and the other a negative coordinate because they straddle the center of the box, the box dimensions must be known at each timestep to get their actual separation. Is there an easier way to handle this type of simulation? I don't think I want to freeze or restrain the groups that need to be fixed, because they must still be allowed flexibility to twist and deform--I only want to roughly fix their centers of mass. Any suggestions would be appreciated, Thanks, Adam -- 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] restraints in PMF (Justin's tutorial)
Actually, this depends on what you are trying to achieve. If you simply want to obtain the standard binding free energy, and somehow you know that the bound state is represented by your umbrella at dist=0 (via a crystal structure, for example), then using additional restraints is common, acceptable, and in many cases desirable. You can do this in analogy to, for example: D. L. Mobley, J. D. Chodera, K. A. Dill. On the use of orientational restraints and symmetry number corrections in alchemical free energy calculations, Journal of Chemical Physics 125:084902, 2006. Selected for the Virtual Journal of Biological Physics Research 12(5), 2006. Since the free energy is a state function, you are free to select any pathway and it is useful to pick one that converges quickly... and those involving rotation of large molecules are sure to converge very slowly. If you do this, then you'll need to account for the restraints in your free energy term in a way similar to that outlined in the paper that I referenced above. The only reason that you would want to stick to the simple PMF is if you are actually interested in the shape of the unrestrained PMF, in which case you can't add additional restraints. Chris. -- original message -- Rebeca García Fandiño wrote: Thanks a lot for your quick answer! I think they are separated enough, however my monomers are cyclic (like discs); I start with a parallel conformation between then, but along the Umbrella simulation, both of them rotate and approach. If I do not use restraints, how could I avoid the rotation? You don't. Why wouldn't two molecules rotate freely in solution when binding or unbinding? It seems like completely natural behavior. Even in simple systems of protein-ligand association, part of the binding energy is the entropic restriction of the ligand into a certain binding-competent pose. Why wouldn't that happen here? Sounds like an artificial notion to me. -Justin I am using the following md_umbrella.mdp: title = Umbrella pulling simulation define = define = ; Run parameters integrator = md dt = 0.002 tinit = 0 nsteps = 50 ; 1 ns nstcomm = 10 ; Output parameters nstxout = 5000nstvout = 5000 nstfout = 5000 nstxtcout = 5000nstenergy = 5000 ; Bond parameters constraint_algorithm= lincs constraints = all-bonds continuation= yes ; Single-range cutoff scheme nstlist = 5 ns_type = grid rlist = 1.4 rcoulomb= 1.4 rvdw= 1.4 ; PME electrostatics parameters coulombtype = PME fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 ewald_rtol = 1e-5 optimize_fft= yes ; Berendsen temperature coupling is on in two groups Tcoupl = Nose-Hoover tc_grps = r_1_r_2 CL3 tau_t = 0.5 0.5 ref_t = 300 300 ; Pressure coupling is on Pcoupl = Parrinello-Rahman pcoupltype = isotropic tau_p = 1.0 compressibility = 4.5e-5 ref_p = 1.0 ; Generate velocities is off gen_vel = no ; Periodic boundary conditions are on in all directions pbc = xyz ; Long-range dispersion correction DispCorr= EnerPres ; Pull code pull= umbrella pull_geometry = distance pull_dim= N N Y pull_start = yes pull_ngroups= 1 pull_group0 = r_1 pull_group1 = r_2 pull_init1 = 0 pull_rate1 = 0.0 pull_k1 = 1000 ; kJ mol^-1 nm^-2 pull_nstxout= 1000 ; every 2 ps pull_nstfout= 1000 ; every 2 ps Thanks a lot again for your help. Best wishes, Rebeca. Date: Wed, 22 Jun 2011 10:53:16 -0400 From: jalemkul at vt.edu To: gmx-users at gromacs.org Subject: Re: [gmx-users] restraints in PMF (Justin's tutorial) Rebeca García Fandiño wrote: Hello, I am trying to obtain the PMF from Umbrella Sampling of the process of separating two monomers of a dimer, following Justin 's tutorial http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/umbrella/index.html I have done the Umbrella Sampling simulations without using any restraints in any of the monomers, however, I can see that they move and gyrate so that although the c.o.m are separated from each other, there are parts of both interacting, in such way that they are not separated as they should be. Would it be correct if I apply restraints to both monomers in all the Umbrella Sampling simulations?. I have seen that in Justin's tutorial they applied restraints to one of the chains, but in my case I think I will need to restrain both of the units. Would that be correct for the PMF calculation? There are no position restraints applied during the umbrella sampling simulations. They are unnecessary. The umbrella potential is itself a restraint to maintain COM separation. If parts of your proteins are
[gmx-users] Timing variability
Dear Users:
Has anybody else looked at simulation speed (ns/day) over the segments
of long runs? I always benchmark and optimize my systems carefully,
but it was only recently that I realized how much variability I am
obtaining over long runs. Perhaps this is specific to my cluster,
which is one reason for my post.
Here is the performance (in ns/day) that I obtain with a single run.
This is representative of what I see with 30 other runs (see list
below). First, the distribution is bimodal, with peaks at around 80
and 90 ns/day. Second, the values go as low as 70 ns/day (and I have
seen as low as 50 ns/day when I look through all 30 run directories
that differ only by the position of the umbrella restraint).
I am using gromacs 4.5.3 and I run it like this:
mpirun mdrun_mpi -deffnm md1 -dlb yes -npme 16 -cpt 30 -maxh 47 -cpi
md1.cpt -cpo md1.cpt -px coord.xvg -pf force.xvg -noappend
I have obtained similar behaviour also with another system.
I have attempted to correlate timing with the node on which the job
runs, the output of the ^Grid: in the .log file, the DD division of
PME and real-space and the value of pme mesh/force, all to no avail. I
have found one correlation whereby the very slow runs also indicate a
high relative load for PME.
I suspect that it is some value that is being determined at run start
time that is affecting my performance, but I am not sure what this
could be. Perhaps the fourier-spacing is leading to different initial
grids?
Thank you (timing information and my .mdp file follows),
Chris
### Output the ns/day obtained by the run segments
$ grep Perf *log|awk '{print $4}'
78.286
82.345
81.573
83.418
92.423
90.863
85.833
91.131
91.820
71.246
76.844
91.805
92.037
85.702
92.251
89.706
88.590
89.381
90.446
81.142
76.365
76.968
76.037
79.286
79.895
79.047
78.273
79.406
78.018
78.645
78.172
80.255
81.032
81.047
77.414
78.414
80.167
79.278
80.892
82.796
81.300
77.392
71.350
73.134
76.519
75.879
80.684
81.076
87.821
90.064
88.409
80.803
88.435
### My .mdp file
constraints = all-bonds
lincs-iter = 1
lincs-order = 6
constraint_algorithm = lincs
integrator = sd
dt = 0.004
tinit = 100
init_step= 0
nsteps = 10
nstcomm = 1
nstxout = 10
nstvout = 10
nstfout = 10
nstxtcout = 25000
nstenergy = 25000
nstlist = 5
nstlog=0 ; reduce log file size
ns_type = grid
rlist = 1
rcoulomb = 1
rvdw = 1
coulombtype = PME
ewald-rtol = 1e-5
optimize_fft = yes
fourierspacing = 0.12
fourier_nx = 0
fourier_ny = 0
fourier_nz = 0
pme_order = 4
tc_grps = System
tau_t = 1.0
ld_seed = -1
ref_t = 300
gen_temp = 300
gen_vel = yes
unconstrained_start = no
gen_seed = -1
Pcoupl = berendsen
pcoupltype = semiisotropic
tau_p = 4 4
compressibility = 4.5e-5 4.5e-5
ref_p = 1.0 1.0
; COM PULLING
pull = umbrella
pull_geometry= position
pull_dim = N N Y
pull_start = no
pull_nstxout = 250
pull_nstfout = 250
pull_ngroups = 1
pull_group0 = POPC
pull_pbcatom0= 338
pull_group1 = KSC
pull_pbcatom1= 0
pull_init1 = 0 0 0.0
pull_rate1 = 0
pull_k1 = 3000.0
pull_vec1= 0 0 0
;;;EOF
--
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] T-A mutation for a dimer protein
Dear Lishan: First, it would be great to see some evidence that you have tried to do this yourself before posting. Your I think makes it a possible waste of time for us to suggest a resolution to a problem that may or may not exist. Second, if indeed it is a problem, perhaps you could use two identical topologies with different molecule names and treat them separately. Chris -- original message -- Hi All, I have a protein dimer and I want to calculate a T to A mutation free energy change using TI method. Since it is a dimer, it is very convenient (and advantageous) to mutate the T in both monomer simultaneously. Gromacs will write out dH/dl sum for the two mutations together, I think. My question is how can I change Gromacs so than dH/dl for T-A mutation in each monomer is written out separately? Best, Lishan -- 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] Programs to add residues
Try Loopy. You can get it to build termini in addition to loops. http://wiki.c2b2.columbia.edu/honiglab_public/index.php/Software:Loopy Nevertheless, I'd suggest simply omitting that part of the protein and capping your new terminus to remove the charge. You will have more difficulties converging the conformation of the unstructured terminus than you may expect. CHris. -- original message -- Hi all, Is there a program that allows the user to add residues to the N and C terminus, without using the electron density? I would like to add a short linker to my protein which doesn't exist in the electron density. Thanks a lot, Sincerely, Zack -- 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] is lincs used with virtual hydrogens?
Thank you Mark, I really appreciate it. -- original message -- On 20/06/2011 1:23 AM, chris.neale at utoronto.ca wrote: Dear Mark: Now I am confused. Your first post indicated that P-LINCS did the angle constraints. But here you indicate that the v-site algorithm does it. No... my second email observed that LINCS warnings can result when v-sites are being used, and I ascribed that effect to the size of the time step. I didn't say v-sites were doing anything with angle constraints. This is probably because my first post was incomplete about the method that I used. Can you please confirm my current understanding? Therefore (A): pdb2gmx -vsite none constraints = h-angles constraint_algorithm = lincs -- P-LINCS constrains all bonds and angles involving hydrogens But when I use (B): pdb2gmx -vsite hydrogen constraints = all-bonds constraint_algorithm = lincs -- Hydrogens are built as V-sites and P-LINCS constrains all bonds that do not involve hydrogen. Of course, water is treated by SETTLE in both of the above. Please note that I am not concerned about getting LINCS warnings. I probably should not have mentioned that. I am simply trying to figure out exactly how to write my methods section. Both the above are correct understandings. I think you misread my second email. I'd intended to put my parenthetical observation about LINCS in the first email, but forgot to (however you seem to have understood my context fine). Sorry for whatever. :) Mark -- 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] is lincs used with virtual hydrogens?
Dear Mark: Now I am confused. Your first post indicated that P-LINCS did the angle constraints. But here you indicate that the v-site algorithm does it. This is probably because my first post was incomplete about the method that I used. Can you please confirm my current understanding? Therefore (A): pdb2gmx -vsite none constraints = h-angles constraint_algorithm = lincs -- P-LINCS constrains all bonds and angles involving hydrogens But when I use (B): pdb2gmx -vsite hydrogen constraints = all-bonds constraint_algorithm = lincs -- Hydrogens are built as V-sites and P-LINCS constrains all bonds that do not involve hydrogen. Of course, water is treated by SETTLE in both of the above. Please note that I am not concerned about getting LINCS warnings. I probably should not have mentioned that. I am simply trying to figure out exactly how to write my methods section. Thank you, Chris. -- original message -- Mark On Sat, Jun 18, 2011 at 4:40 PM, chris.neale at utoronto.ca mailto:chris.neale at utoronto.ca chris.neale at utoronto.ca mailto:chris.neale at utoronto.ca wrote: Thank you Roland. I did use: constraints = all-bonds lincs-iter = 1 lincs-order = 6 constraint_algorithm = lincs From looking at the manual, I figured that angle and bond constraints would all be done by LINCS if I had done (A): pdb2gmx -vsite none constraints = h-angles (a combination that I have never tried) But when I use (B): pdb2gmx -vsite hydrogen constraints = all-bonds It seems possible to me that LINCS is not used but instead the position of the atom is simply built from a mathematical function. Perhaps this all stems from my lack of thorough understanding of LINCS, but it seems to me that there need be no iteration to simply place an atoms based on virtual_sites3 (which are constructed by pdb2gms -hydrogen) For now, I'll simply add a line to state that I built virtual sites for hydrogen atoms to make it clear, but I'd still like to understand the difference between options A and B, above, if you have some time. Thank you again, Chris. On Sat, Jun 18, 2011 at 4:13 PM, chris.neale at utoronto.ca http://utoronto.ca chris.neale at utoronto.ca http://utoronto.ca wrote: Dear Users: If I create the topology of a peptide like this: pdb2gmx -f protein.gro -vsite hydrogens And then simulate it in vacuum, is lincs used at all? I believe that it is, as if I use a timestep that is too large then I get LINCS warnings about angles rotating more than 30 degrees, but that warning message could possibly have been written with the assumption that I used LINCS and not virtual hydrogens. Probably. To make sure check the constraint-algorithm selected in your mdp. BTW: If you want to use large timecheck you should normally use constraints=all-bonds and lincs-order=6. Finally, is there a method that needs to be named or cited in relation to the fact that the angles are now constrained? Is that also done with P-LINCS? This is also done with P-LINCS. Not sure whether one should sign something regarding the construction/usage of v-sites. Roland Thank you, Chris. -- 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] is lincs used with virtual hydrogens?
Dear Users: If I create the topology of a peptide like this: pdb2gmx -f protein.gro -vsite hydrogens And then simulate it in vacuum, is lincs used at all? I believe that it is, as if I use a timestep that is too large then I get LINCS warnings about angles rotating more than 30 degrees, but that warning message could possibly have been written with the assumption that I used LINCS and not virtual hydrogens. Finally, is there a method that needs to be named or cited in relation to the fact that the angles are now constrained? Is that also done with P-LINCS? Thank you, Chris. -- 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] is lincs used with virtual hydrogens?
Thank you Roland. I did use: constraints = all-bonds lincs-iter = 1 lincs-order = 6 constraint_algorithm = lincs From looking at the manual, I figured that angle and bond constraints would all be done by LINCS if I had done (A): pdb2gmx -vsite none constraints = h-angles (a combination that I have never tried) But when I use (B): pdb2gmx -vsite hydrogen constraints = all-bonds It seems possible to me that LINCS is not used but instead the position of the atom is simply built from a mathematical function. Perhaps this all stems from my lack of thorough understanding of LINCS, but it seems to me that there need be no iteration to simply place an atoms based on virtual_sites3 (which are constructed by pdb2gms -hydrogen) For now, I'll simply add a line to state that I built virtual sites for hydrogen atoms to make it clear, but I'd still like to understand the difference between options A and B, above, if you have some time. Thank you again, Chris. On Sat, Jun 18, 2011 at 4:13 PM, chris.neale at utoronto.ca chris.neale at utoronto.ca wrote: Dear Users: If I create the topology of a peptide like this: pdb2gmx -f protein.gro -vsite hydrogens And then simulate it in vacuum, is lincs used at all? I believe that it is, as if I use a timestep that is too large then I get LINCS warnings about angles rotating more than 30 degrees, but that warning message could possibly have been written with the assumption that I used LINCS and not virtual hydrogens. Probably. To make sure check the constraint-algorithm selected in your mdp. BTW: If you want to use large timecheck you should normally use constraints=all-bonds and lincs-order=6. Finally, is there a method that needs to be named or cited in relation to the fact that the angles are now constrained? Is that also done with P-LINCS? This is also done with P-LINCS. Not sure whether one should sign something regarding the construction/usage of v-sites. Roland Thank you, Chris. -- 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] gmx4.5.4 genion problem: No line with moleculetype 'SOL' found
For the quick fix: 1. run genion on your topology that does work. Look at this to see the format of the ion atom and residue names 2. Pick a few waters in the structure containing the ligand and replace the OW by the ion and remove the hydrogens, then fix the number of atoms on the second line of the file and run editconf on the file. ** But I wouldn't be satisfied with that if I were you. there is some problem here that you should uncover now before you move farther. Show us your commands exactly (copy and paste) and perhaps we can help more. As a side note: your topology that you presented is going to cause you problems if you try to put position restraints on the water, as it is incorrectly ordered: ; Include water topology #include gromos43a1.ff/spc.itp ; Include ligand topoloty #include drg.itp #ifdef POSRES_WATER ; Position restraint for each water oxygen [ position_restraints ] ; i funct fcxfcyfcz 11 1000 1000 1000 #endif -- original message -- Thank you for your suggestions. I have checked the topology file with vim, and it looks pefectly, also, the only problem happens when I use genion. One thing that might be possible is that I use the double precision version of gromacs, because when I solve it in water, the written of topology file looks weired(like I showed). What do you mean by work-around? You mean just manually change some water molecules in topology file and coordinate file into ions? Will the program recognize something like CL, NA? Or can I just add ions to the system and solve it in water? Thank you very much Best Wishes -- 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] g_dist
Dear Nilesh: You don't seem to have made any progress since Mark gave you some reading hints. Your best chance to get a useful response follows the general form: 1. I want to do A. 2. I tried method B, here are the exact commands that I used (copy and paste) 3. With that method, I obtain C 4. But I am confused about why C does not look like D. 5. What is the problem? Basically, it helps to show that you are doing your own work too. For example, try a command on a single frame and then look at that frame with VMD and count for yourself and verify the answer that way. Chris. -- original message -- I have a system with 128 emi (cations) and 128 Cl (anions). I want to calculate how many CL atoms are in cutoff distance relative to hydorgen aotm of cation. I considered all CL atoms are distinguishable. Basically I want to calcualte the distance between each CL atom and correponing hydorgen and registered those CL atoms which are whthin cutoff. How can I do ? I am using Gromacs 4.0.7 version. Thanks Nilesh -- 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] Can't unfold the protein
Dear Hsin-Lin: The first thing that comes to mind is using very high temperatures with NVT. Obtaining the desired denatured state is a complex challenge, but you might try this paper: Phys. Rev. Lett. 93, 238105 (2004) Reversible Temperature and Pressure Denaturation of a Protein Fragment: A Replica Exchange Molecular Dynamics Simulation Study Still, your post seems to be about simply being unable to unfold at all. I have unfolded a very stable beta-barrel by either simulating at 3000 K or using the pull code to pull the termini apart. This was not a real study -- I only used this to make a reverse folding movie for non-scientists -- but I can verify that it is possible to unfold even very stable proteins by these methods. Chris. -- original message -- Hi, I have question about unfolding. I use three different ways respectively but all failed. The three way is, 1. 600K in 20ns, then the system explode. 2. 400K in 10ns, the protein is very stable and the value is almost the same in radius of gyration. 3. heating up from 300K to 400K in 2ns and the other 2ns for 400K only, then the protein is still stable and the value is almost the same in radius of gyration. Below is the mdp file of the heating. What's wrong with my system? -- 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] Can't unfold the protein
Dear Hsin-Lin: I am no expert in this area. I am just saying that if you get densities that are way too low in NPT, then you might alleviate this problem with NVT. In fact, I would personally do this in vacuum without pbc and use the sd integrator. Then you can sample extended conformations. I'd consider this a method to generate unfolded conformations that probably have no relation to the true temperature denatured state. But let's be honest, you're never going to come close to Boltzmann sampling of an unfolded state in 300 ns so why bother with the water? Your link ( http://manual.gromacs.org/online/protunf.html ) is only concerned with analysis, so it is irrelevant. If you are determined to use NPT then I suppose that you could use the Berendsen barostat (more stable in simulations but you get the wrong ensemble) with a very low compressibility. That's about all I can say, perhaps somebody else will comment. Good luck, Chris. -- original message -- Dear Chris, Thank you for your reply. My protein is very stable. I simulated it in 300ns in 300K before but there was almost no change. That's why I want to do denaturation now. I want to start a new simulation on the protein which is unfolded. And I'll read the paper you told me. Thank you. But I don't understand what you recommend me to do. You use 3000K on your protein and the protein unfolded. I use 600K on my protein but the system explode(box become very big and protein locate at corner). I saw the second way on this web, http://manual.gromacs.org/online/protunf.html But it also useless to me. And in my mdp I use thermostat and barostat, which means my system is NPT. Does anything wrong in my methods? Sincerely yours, Hsin-Lin -- 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] g_hbond
From g_hbond -h OH and NH groups are regarded as donors, O is an acceptor always, N is an acceptor by default, but this can be switched using -nitacc. Dummy hydrogen atoms are assumed to be connected to the first preceding non-hydrogen atom. The two groups are to find hbonds between group A and group B. If you want to find hbonds within group A, then you specify an index group that contains group A two times. I think that you want: [ group ] 100 101 200 Then give group 2x as the identical index group. PS, there were problems with g_hbond a few months ago reported on list, I suggest that you check to ensure that these problems were fixed. Chris. -- original message -- Hi, To help clarify my previous g_hbond question, I have these two groups in my index file [ a_100_a_101] # 100 is N, #101 is H 100 101 [ a_200 ]# 200 is O 200 when I ran g_hbond, I got Found 1 donors and 2 acceptors. The part that I don't quite understand is why it says 2 acceptors. Did I set up the groups correctly? Again, hope this clarifies my question. Thanks for your insight. Simon -- 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] centering a bilayer at a cartesian coordinate
Dear users:
I am trying to process a trajectory so that a group of molecules has
their center of mass at a constant position in Cartesian space. I have
not been able to figure out how to do this.
The reason that I would like this is that I have conducted umbrella
sampling of a solute along the normal to a lipid bilayer and I would
like to create a movie of a false-trajectory of the immersion profile.
Such processing is also important for generating SDFs and for use with
g_density.
I expect that the problems that I am having are due to volume
fluctuations in the NPT ensemble.
The best idea that I have tried was (gromacs 4.5.3):
trjconv -center -box 10 10 15
where the original dimensions are much less that 10 10 15
But still, when I run g_traj on the processed .gro file I do not get
similar values for the COM.
In more detail:
echo 0 | trjconv -s md1.tpr -f md1_50ps_md1.part0001.xtc -o tmp.xtc -e
10 -pbc mol
# select POPC for centering and System for output and use a .tpr with pbc=no
echo -e 13\n0\n | trjconv -s empty.tpr -f tmp.xtc -box 10 10 15
-center -pbc none -o tmp2.xtc
# select POPC for coordinate output
echo 13 | g_traj -s empty.tpr -f tmp2.xtc -ox -com -nopbc
99550 4.8424 5.00025 7.49519
99600 4.79691 4.94304 7.39868
99650 4.75789 4.85189 7.41581
99700 4.74974 4.80423 7.48257
99750 5.07858 5.04788 7.51654
99800 4.99091 4.78976 7.47778
99850 5.11406 4.84656 7.48769
99900 5.09739 5.12395 7.47889
99950 5.0411 4.88659 7.52737
10 4.91271 5.07101 7.50397
For Z, the value is 7.48037+/-0.072562 (MIN=7.3 and MAX=7.8)
If I look at the frames from the min and max z values, they do clearly
differ in the placement of the bilayer along z.
###
To simplify things, I redid this while only using the phosphorous
atoms (1 per lipid) in the centering and coordinate output.
Here, Z=7.49323+/-0.0658727 (MIN=7.3 and MAX=7.7)
To verify this and since all the atoms have the same mass, I checked
it with awk:
echo 0 | trjconv -s empty.tpr -f tmp2.xtc -dump 47550 -b 4 -o look_low.gro
echo 0 | trjconv -s empty.tpr -f tmp2.xtc -dump 80750 -b 8 -o
look_high.gro
cat look_low.gro|grep P8|awk '{s+=$NF;n++} END{print s/n}'
7.39197
cat look_high.gro|grep P8|awk '{s+=$NF;n++} END{print s/n}'
7.62207
###
Finally, I rechecked with a single atom and here it worked exactly as
expected. The atoms is placed at 5 5 7.5 in every single frame.
###
Then I checked using 2 phosphorous atoms from the same (or,
separately, different) leaflets, and the result was the same:
$ tail coord.xvg
99550 5 5 7.5
99600 5 5.0005 7.5
99650 5.0005 5.0005 7.5005
99700 5 5.0005 7.5
99750 5.0005 5 7.5005
99800 5 5.0005 7.5005
99850 5 5 7.5005
99900 5.0005 5 7.5005
99950 5.0005 5 7.4995
10 5 5 7.5
Then using 3 phosphorous atoms from (containing atoms from both
leaflets), and now there is more deviation:
$ tail coord.xvg
99550 5.48833 5.30233 6.972
99600 5.48733 5.29233 6.97167
99650 5.50033 5.29567 6.93467
99700 5.52133 5.27433 6.93967
99750 5.489 5.253 6.99133
99800 5.481 5.319 6.91533
99850 5.47867 5.309 6.827
99900 5.48733 5.292 6.83133
99950 5.52767 5.276.86867
10 5.518 5.319 6.89167
(MIN=6.7 and MAX=7.1)
I originally had problems with pbc during trjconv in which the -center
option would place the bilayer at the peripheries of the unit cell in
some frames (COM should still be at the center). This is what lead me
to the route that I have used above (a .tpr that was generated with
pbc=no and trjconv -pbc none -center). Nevertheless, I think that
perhaps trjconv is using pbc information even in this case where ti
should not.
##
If anybody can see what I am doing wrong or has a method to center a
bilayer at a cartesian coordinate, then I am very interested.
Thank you for your time,
Chris.
--
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] centering a bilayer at a cartesian coordinate
I solved it. It was my own error to assume that trjconv -center
centers the COM. Actually, it centers the value of (min+max)/2 in each
dimension.
I can modify trjconv locally to do what I want.
Thank you,
Chris.
static void center_x(int ecenter,rvec x[],matrix box,int n,int
nc,atom_id ci[])
{
int i,m,ai;
rvec cmin,cmax,box_center,dx;
if (nc 0) {
copy_rvec(x[ci[0]],cmin);
copy_rvec(x[ci[0]],cmax);
for(i=0; inc; i++) {
ai=ci[i];
for(m=0; mDIM; m++) {
if (x[ai][m] cmin[m])
cmin[m]=x[ai][m];
else if (x[ai][m] cmax[m])
cmax[m]=x[ai][m];
}
}
calc_box_center(ecenter,box,box_center);
for(m=0; mDIM; m++)
dx[m] = box_center[m]-(cmin[m]+cmax[m])*0.5;
for(i=0; in; i++)
rvec_inc(x[i],dx);
}
}
Quoting chris.ne...@utoronto.ca:
Dear users:
I am trying to process a trajectory so that a group of molecules has
their center of mass at a constant position in Cartesian space. I have
not been able to figure out how to do this.
The reason that I would like this is that I have conducted umbrella
sampling of a solute along the normal to a lipid bilayer and I would
like to create a movie of a false-trajectory of the immersion profile.
Such processing is also important for generating SDFs and for use with
g_density.
I expect that the problems that I am having are due to volume
fluctuations in the NPT ensemble.
The best idea that I have tried was (gromacs 4.5.3):
trjconv -center -box 10 10 15
where the original dimensions are much less that 10 10 15
But still, when I run g_traj on the processed .gro file I do not get
similar values for the COM.
In more detail:
echo 0 | trjconv -s md1.tpr -f md1_50ps_md1.part0001.xtc -o tmp.xtc -e
10 -pbc mol
# select POPC for centering and System for output and use a .tpr with pbc=no
echo -e 13\n0\n | trjconv -s empty.tpr -f tmp.xtc -box 10 10 15
-center -pbc none -o tmp2.xtc
# select POPC for coordinate output
echo 13 | g_traj -s empty.tpr -f tmp2.xtc -ox -com -nopbc
99550 4.8424 5.00025 7.49519
99600 4.79691 4.94304 7.39868
99650 4.75789 4.85189 7.41581
99700 4.74974 4.80423 7.48257
99750 5.07858 5.04788 7.51654
99800 4.99091 4.78976 7.47778
99850 5.11406 4.84656 7.48769
99900 5.09739 5.12395 7.47889
99950 5.0411 4.88659 7.52737
104.91271 5.07101 7.50397
For Z, the value is 7.48037+/-0.072562 (MIN=7.3 and MAX=7.8)
If I look at the frames from the min and max z values, they do clearly
differ in the placement of the bilayer along z.
###
To simplify things, I redid this while only using the phosphorous
atoms (1 per lipid) in the centering and coordinate output.
Here, Z=7.49323+/-0.0658727 (MIN=7.3 and MAX=7.7)
To verify this and since all the atoms have the same mass, I checked
it with awk:
echo 0 | trjconv -s empty.tpr -f tmp2.xtc -dump 47550 -b 4 -o
look_low.gro
echo 0 | trjconv -s empty.tpr -f tmp2.xtc -dump 80750 -b 8 -o
look_high.gro
cat look_low.gro|grep P8|awk '{s+=$NF;n++} END{print s/n}'
7.39197
cat look_high.gro|grep P8|awk '{s+=$NF;n++} END{print s/n}'
7.62207
###
Finally, I rechecked with a single atom and here it worked exactly as
expected. The atoms is placed at 5 5 7.5 in every single frame.
###
Then I checked using 2 phosphorous atoms from the same (or,
separately, different) leaflets, and the result was the same:
$ tail coord.xvg
99550 5 5 7.5
99600 5 5.0005 7.5
99650 5.0005 5.0005 7.5005
99700 5 5.0005 7.5
99750 5.0005 5 7.5005
99800 5 5.0005 7.5005
99850 5 5 7.5005
99900 5.0005 5 7.5005
99950 5.0005 5 7.4995
105 5 7.5
Then using 3 phosphorous atoms from (containing atoms from both
leaflets), and now there is more deviation:
$ tail coord.xvg
99550 5.48833 5.30233 6.972
99600 5.48733 5.29233 6.97167
99650 5.50033 5.29567 6.93467
99700 5.52133 5.27433 6.93967
99750 5.489 5.253 6.99133
99800 5.481 5.319 6.91533
99850 5.47867 5.309 6.827
99900 5.48733 5.292 6.83133
99950 5.52767 5.276.86867
105.518 5.319 6.89167
(MIN=6.7 and MAX=7.1)
I originally had problems with pbc during trjconv in which the -center
option would place the bilayer at the peripheries of the unit cell in
some frames (COM should still be at the center). This is what lead me
to the route that I have used above (a .tpr that was generated with
pbc=no and trjconv -pbc none -center). Nevertheless, I think that
perhaps trjconv is using pbc information even in this case where ti
should not.
##
If anybody can see what I am doing wrong or has a method to center a
bilayer at a cartesian coordinate, then I am very interested.
Thank you for your time,
Chris.
--
gmx-users
[gmx-users] g_density
Dear Matthias: Did you run trjconv -center -pbc mol at some point before running g_density? With pbc, there are multiple ways to put the center of a system at the center of the box. This can lead to having the water slab in the center of the unit cell and the bilayer on the top and bottom z. Without seeing the commands and the output, it is hard to say any more. Chris. -- original message -- Hi, I have aligned a trajectory to a reference structure and have now run g_density on both the original and the aligned trajectory in order to find the bilayer headgroup position and the bilayer thickness. g_density works fine with the original trajectory but there seems to be a bug when running it to the aligned structure (Alignment done using trjconv -fit progressive). The density is calculated along the z direction (-d Z). The density for positive z is calculated correctly, but the density for negative z is shifted about 15nm to to the top end of the box size, as if the density distribution along the z-axis were discontinuous. Other then stitching the data together by hands, are there any obvious better solutions? Thanks, Matthias -- 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] Question about umbrella sampling through a nanotube
Unfortunately, I don't think that there is any way to use this data (and only this data) to derive the PMF along z. The way that you did your US, the PMF along z is convoluted with xy motion. You can get the PMF along the reaction coordinate that you actually used (XYZ) using standard WHAM and try to interpret it, I suppose. In your special case, there is probably a mathematical conversion between the XYZ PMF and the Z-only PMF due to the symmetry of the system. What one would desire is two separate restraints, one operating in Z and one operating in XY. That is, unfortunately, not possible in the standard version of gromacs. -- original message -- Dear GMXers, I'm trying to compute the PMF of a molecule along the centerline of a nanotube (the axial direction of the nanotube is parallel to the z axis). The nanotube is used as the reference group and the molecule as the pulling group. pull_geometry = position pull_dim= Y Y Y pull_start = no pull_init1 = 0 0 1.2 pull_rate1 = 0.00 pull_k1 = 1000 pull_nstxout= 1000 pull_nstfout= 1000 pull_dim=Y Y Y is used to make sure that the molecule is always close to the centerline of the nanotube. From the output of pullx.xvg, the COM positions in all three directions (i.e. x,y,z) are included. How to get the PMF only along z axis from here? Is it OK if I just take the corresponding z column and do g_wham directly? Thank you. Best, Yanbin -- next part -- -- 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] Gradually increasing the force constant during restrained MD runs
Dear Aditi:
you can do this by scripting gromacs to run in short segments,
modifying the force constant for each run. e.g.:
#!/bin/bash
for((i=1;i=100; i++)); do
let j=i-1
grompp -f ${i}.mdp -p ${i}.top -f ${j}.gro -o ${i}.tpr
mdrun -deffnm ${i}
done
You only need to use some more scripting to create the ${i}.top and
${i}.mdp files, although I think that orientational restraints are
going to be part of the topology so the .mdp file can even remain
constant.
You can then trjcat -settime once you are done.
-- original message --
I am running an MD simulation on RNA with orientation restraints.
Instead of using a single force constant for the whole trajectory, I
want to start with a very low force constant and increase it gradually
and regularly up to a final suitably high value and then continue the
simulation at this final value of force constant.
Can you suggest if there is a way to do so in GROMACS?
Thanks!
Aditi
--
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] gromacs 2.1
Dear users: Does anybody have the source code for gromacs version 2.1? I would like to check the original source of g_order, but only versions 3.0 and above are available on the gromacs website. Thank you, Chris. -- 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] Calculation of unsaturated deuterium order parameters for POPC
I have patched g_order version 4.5.4 to compute the order parameters
in a different way. Despite trying, I don't understand the original
g_order implementation. In the version that I am providing, the
hydrogens are explicitly built and the coordinates of these hydrogens
are output to stdout.
Compared to the original g_order: there are large differences for the
2 carbons in the double bond and all single bonded carbons are the
same within what appears to be rounding error.
Compared to the VMD script that Thomas Piggot mentioned, the results
differ more, but only because that VMD script is using the real
hydrogen positions and not some with idealized geometry. There is a
larger difference for the double-bonded carbons, but if I use the
idealized hydrogen positions (output by my new version of g_order)
with the VMD .tcl script then I get the same values.
notes on the patch:
1. it is not complete for all options. Use it only with g_order -od
and with or without -unsat.
2. There are extraneous print statement that will provide you with the
coordinates of the atoms that are built. This is for debugging. Either
run with g_order /dev/null or comment out the 3 fprintf lines near
the comments containing the word LOOK in capitals.
3. I did my best, but it could always be in error. Please use with caution.
To apply the patch (version 4.5.4)
cd src/tools
patch g_order.patch
If it will not compile, then the spacing/hard-returns of the patch
must have been mangled in the process of posting. Contact me off-list
and I will provide you with a copy.
Here is the patch:
--- gmx_order.c 2011-03-04 06:10:44.0 -0500
+++ gmx_order.c 2011-06-09 17:57:11.586591000 -0400
@@ -379,6 +379,13 @@
real arcdist;
gmx_rmpbc_t gpbc=NULL;
+// BEGIN CN MOD
+// variables for modification by CN
+rvec ab,ac,bc,e,ce,eg,ef,g,f,cg,cf;
+float dab,def,deg,gdot,fdot,edot;
+// END CN MOD
+
+
/* PBC added for center-of-mass vector*/
/* Initiate the pbc structure */
memset(pbc,0,sizeof(pbc));
@@ -501,98 +508,66 @@
direction[0],direction[1],direction[2]);*/
}
- if (bUnsat) {
- /* Using convention for unsaturated carbons */
- /* first get Sz, the vector from Cn to Cn+1 */
- rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], dist);
- length = norm(dist);
- check_length(length, a[index[i]+j], a[index[i+1]+j]);
- svmul(1/length, dist, Sz);
-
- /* this is actually the cosine of the angle between the double bond
-and axis, because Sz is normalized and the two other components of
-the axis on the bilayer are zero */
- if (use_unitvector)
- {
- sdbangle += acos(iprod(direction,Sz)); /*this can probably be
optimized*/
- }
- else
- sdbangle += acos(Sz[axis]);
- } else {
- /* get vector dist(Cn-1,Cn+1) for tail atoms */
- rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i-1]+j]], dist);
- length = norm(dist); /* determine distance between two atoms */
- check_length(length, a[index[i-1]+j], a[index[i+1]+j]);
-
- svmul(1/length, dist, Sz);
- /* Sz is now the molecular axis Sz, normalized and all that */
- }
-
- /* now get Sx. Sx is normal to the plane of Cn-1, Cn and Cn+1 so
- we can use the outer product of Cn-1-Cn and Cn+1-Cn, I hope */
- rvec_sub(x1[a[index[i+1]+j]], x1[a[index[i]+j]], tmp1);
- rvec_sub(x1[a[index[i-1]+j]], x1[a[index[i]+j]], tmp2);
- cprod(tmp1, tmp2, Sx);
- svmul(1/norm(Sx), Sx, Sx);
-
- /* now we can get Sy from the outer product of Sx and Sz */
- cprod(Sz, Sx, Sy);
- svmul(1/norm(Sy), Sy, Sy);
-
- /* the square of cosine of the angle between dist and the axis.
- Using the innerproduct, but two of the three elements are zero
- Determine the sum of the orderparameter of all atoms in group
- */
- if (use_unitvector)
- {
- cossum[XX] = sqr(iprod(Sx,direction)); /* this is allowed, since Sa
is normalized */
- cossum[YY] = sqr(iprod(Sy,direction));
- cossum[ZZ] = sqr(iprod(Sz,direction));
- }
- else
- {
- cossum[XX] = sqr(Sx[axis]); /* this is allowed, since Sa is normalized
*/
- cossum[YY] = sqr(Sy[axis]);
- cossum[ZZ] = sqr(Sz[axis]);
+// BEGIN CN MOD
+ if (bUnsat) {
+//exactly as for the case without unsat, but here we select the
ce vector to represent the hydrogen
+// place the 2 imaginary hydrogen atoms by enscribing a
tetrahedron in a cube
+// a,b,c are the 3 backbone atoms, where a and b are in two cube
corners
[gmx-users] Calculation of unsaturated deuterium order parameters for POPC
Dear Thomas: I agree with your conclusions about g_order. My usage of the Berger POPC lipids in combination with the g_order command (as per the instructions in http://www.gromacs.org/Documentation/Gromacs_Utilities/g_order ) yields order parameters of 0.193 and 0.051 for the first and second carbons in the double bond at 300K. This is similar to but, as expected, slightly larger than the values described in http://www.sciencedirect.com/science/article/pii/S0014579305014365 where they simulated at 310K and also used g_order. I suspect that one could find tens of publications that use g_order and have published similar values. The Poger paper (On the Validation of Molecular Dynamics Simulations of Saturated and cis-Monounsaturated Phosphatidylcholine Lipid Bilayers: A Comparison with Experiment) gives clearly different values for the first carbon in the double bond than is reported by g_order. One probable reason that this has gone unnoticed so long is that there appears to be no experimental oleoyl order parameters -- as far as I know -- if somebody finds them, please post them. I took a look at the source and it is clearly incorrect for the i+1 unsaturated partner in a double bond over i-i+1. I have developed a new version of g_order to fix this, but I'd like to test it first. If we could agree on a single .xtc file to analyze and compare our values then that would be great. I can provide you with one of mine (Berger lipids) or you can provide me with one of yours (off-list). Thanks, Chris. -- 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] Why does the -append option exist?
Dear Dimitar:
to me, your posts have indeed appeared very aggressive. So much so
that there are probably many people who can help who have decided not
to bother due to what they perceive your tone to be.
On to the problem at hand: I agree with the previous suggestion that
you may have somehow had two runs going in the same directory and
competing to write to the same file.
Can you test this by submitting one script, then waiting for an hour,
then submitting the same script again in the same directory? If you
can reproduce the problem with this test then it would suggest that
this is indeed one possible source of the problem, however it arrived.
I probably won't respond on this topic again, as I don't want to find
myself involved in an argument on the list. I simply wanted to point
out that there has been at least one valid suggestion that you are in
a good position to test.
Very generally, when you mail the list, people are going to try to
suggest possible sources of the problem and then, ideally, you can
test those and report back.
Chris.
Here is an example:
[dpachov]$ ll -rth run1* \#run1*
-rw-rw-r-- 1 dpachov dpachov 11K May 2 02:59 run1.po.mdp
-rw-rw-r-- 1 dpachov dpachov 4.6K May 2 02:59 run1.grompp.out
-rw-rw-r-- 1 dpachov dpachov 3.5M May 13 19:09 run1.gro
-rw-rw-r-- 1 dpachov dpachov 2.3M May 14 00:40 run1.tpr
-rw-rw-r-- 1 dpachov dpachov 2.3M May 14 00:40 run1-i.tpr
-rw-rw-r-- 1 dpachov dpachov0 May 29 21:53 run1.trr
-rw-rw-r-- 1 dpachov dpachov 1.2M May 31 10:45 run1.cpt
-rw-rw-r-- 1 dpachov dpachov 1.2M May 31 10:45 run1_prev.cpt
-rw-rw-r-- 1 dpachov dpachov0 Jun 3 14:03 run1.xtc
-rw-rw-r-- 1 dpachov dpachov0 Jun 3 14:03 run1.edr
-rw-rw-r-- 1 dpachov dpachov 15M Jun 3 17:03 run1.log
Submitted by:
ii=1
ifmpi=mpirun -np $NSLOTS
if [ ! -f run${ii}-i.tpr ];then
cp run${ii}.tpr run${ii}-i.tpr
tpbconv -s run${ii}-i.tpr -until 20 -o run${ii}.tpr
fi
k=`ls md-${ii}*.out | wc -l`
outfile=md-${ii}-$k.out
if [[ -f run${ii}.cpt ]]; then
$ifmpi `which mdrun` -s run${ii}.tpr -cpi run${ii}.cpt -v -deffnm
run${ii} -npme 0 $outfile 21
fi
=
This script is not using mdrun -append.
-append is the default, it doesn't need to be explicitly listed.
Your original post suggested the use of -append was a problem. Why aren't
we seeing a script with mdrun -append? Also, please provide the full script
- it looks like there might be a loop around your tpbconv-then-mdrun
fragment.
There is no loop; this is a job script with PBS directives. The header of it
looks like:
===
#!/bin/bash
#$ -S /bin/bash
#$ -pe mpich 8
#$ -ckpt reloc
#$ -l mem_total=6G
===
as usual submitted by:
qsub -N myjob.q
Note that a useful trouble-shooting technique can be to construct your
command line in a shell variable, echo it to stdout (redirected as suitable)
and then execute the contents of the variable. Now, nobody has to parse a
shell script to know what command line generated what output, and it can be
co-located with the command's stdout.
I somewhat understand your point, but could give an example if you think it
is really necessary?
snip
Writing checkpoint, step 51879590 at Tue May 31 10:45:22 2011
Energies (kJ/mol)
U-BProper Dih. Improper Dih. CMAP Dih. LJ-14
8.33208e+034.72300e+035.31983e+02 -1.21532e+032.89586e+03
Coulomb-14LJ (SR) Disper. corr. Coulomb (SR) Coul. recip.
3.00900e+049.31785e+04 -3.87790e+03 -7.40841e+05
-8.36838e+04
PotentialKinetic En. Total EnergyTemperature Pres. DC (bar)
-6.89867e+051.28721e+05 -5.61146e+052.99472e+02 -1.24229e+02
Pressure (bar) Constr. rmsd
-1.03491e+022.99840e-05
So the -append restart looks like it did fine here.
Last output files from restarts:
[dpachov]$ ll -rth md-1-*out | tail -10
-rw-rw-r-- 1 dpachov dpachov 6.1K Jun 3 16:40 md-1-2428.out
-rw-rw-r-- 1 dpachov dpachov 6.2K Jun 3 16:44 md-1-2429.out
-rw-rw-r-- 1 dpachov dpachov 5.9K Jun 3 16:46 md-1-2430.out
-rw-rw-r-- 1 dpachov dpachov 5.9K Jun 3 16:48 md-1-2431.out
-rw-rw-r-- 1 dpachov dpachov 6.1K Jun 3 16:50 md-1-2432.out
-rw-rw-r-- 1 dpachov dpachov0 Jun 3 16:52 md-1-2433.out
-rw-rw-r-- 1 dpachov dpachov 6.2K Jun 3 16:55 md-1-2434.out
-rw-rw-r-- 1 dpachov dpachov 6.2K Jun 3 16:58 md-1-2435.out
-rw-rw-r-- 1 dpachov dpachov 5.9K Jun 3 17:03 md-1-2436.out
*-rw-rw-r-- 1 dpachov dpachov 5.8K Jun 3 17:04 md-1-2437.out*
+ around the time when the run1.xtc file seems to have been saved:
[dpachov]$ ll -rth md-1-23[5-6][0-9]*out
-rw-rw-r-- 1 dpachov dpachov 6.2K Jun
[gmx-users] problem during energy minimization
1. you're not doing energy minimization as your title suggests, please use a better title: integrator = md 2. generally, when you have a problem it is useful to try to simplify the system (e.g. get rid of those tables) and try to reproduce the problem with the simplest system possible before posting. 3. I know that you commented it out, but the following command looks like a typo or a misunderstanding... perhaps you meant -DFLEXIBLE ? ;define= -DEFLEXIBLE 4. Most relevant to you: This is a little ambitious, don't you think? dt = 0.02 ; ps ! Chris. -- original message -- Dear Gmx users, I am trying to do minimization of my system .i have no problem wehen i grompp it but when i do the mdrun its giving me some segmentation error.I had attached the output of grompp and mdrun below.Any suggestions please.Thanks in advance *OUTPUT OF GROMPP* Ignoring obsolete mdp entry 'title' Ignoring obsolete mdp entry 'cpp' Back Off! I just backed up mdout.mdp to ./#mdout.mdp.2# Generated 332520 of the 332520 non-bonded parameter combinations Generating 1-4 interactions: fudge = 0.5 Generated 332520 of the 332520 1-4 parameter combinations Excluding 3 bonded neighbours molecule type 'Ion' Excluding 2 bonded neighbours molecule type 'SOL' Analysing residue names: There are: 2Ion residues There are: 875 Water 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 5253.00 Largest charge group radii for Van der Waals: 0.039, 0.039 nm Largest charge group radii for Coulomb: 0.085, 0.085 nm This run will generate roughly 7 Mb of data Back Off! I just backed up em.tpr to ./#em.tpr.2# *OUTPUT OF MDRUN* Reading file em.tpr, VERSION 4.5.3 (single precision) Starting 2 threads WARNING: For the 1498 non-zero entries for table 2 in table_Na_Cl.xvg the forces deviate on average -2147483648% from minus the numerical derivative of the potential WARNING: For the 1498 non-zero entries for table 2 in table_Na_Cl.xvg the forces deviate on average -2147483648% from minus the numerical derivative of the potential Making 1D domain decomposition 2 x 1 x 1 starting mdrun 'NA SODIUM ION in water' 1 steps,200.0 ps. step 0: Water molecule starting at atom 2553 can not be settled. Check for bad contacts and/or reduce the timestep if appropriate. step 0: Water molecule starting at atom 2595 can not be settled. Check for bad contacts and/or reduce the timestep if appropriate. Wrote pdb files with previous and current coordinates Wrote pdb files with previous and current coordinates Segmentation fault *EM.MDP* title = nacl cpp = /usr/bin/cpp ; the c preprocessor ;define= -DEFLEXIBLE integrator = md dt = 0.02 ; ps ! nsteps = 1 nstlist = -1 coulombtype = user energygrps = Na Cl Sol energygrp_table = Na Cl vdwtype= user ns_type = grid rlist = 1.4 rcoulomb = 1.0 rvdw = 1.0 fourierspacing = 0.12 fourier_nx = 0 fourier_ny = 0 fourier_nz = 0 pme_order = 4 pbc=xyz; ; Energy minimizing stuff emtol = 1000.0 emstep = 0.01 regards, sree -- next part -- -- 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] How to get angle distribution between two tyrosine stacking residues
1. Determine some mathematical calculation that you want to apply to some atoms. 2. make a .ndx group that includes all of those atoms 3. Run g_traj -ox to output the coordinates of those atoms 4. Apply your mathematical calculation using awk. -- original message -- Can anyone tell how to calculate the normal-normal angle between two stacking tyrosine residues as a function of time. -- 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] tabulated potential energy=nan for r=0 nm and charged atoms
Dear users: with assistance from Berk on the developers list, I am posting a work-around to this issue. the nonbonded interactions calculate the nonbonded distance, r, using the gmx_invsqrt function. Thus, the optimized kernel from 4.5.4 and 4.0.7 both give nan when two charges occupy the space, even when using tables that specify the force and potential energy to be zero. I had expected that I would need to modify the generic kernel source, but that turned out to be unnecessary -- just using it was enough. For some reason, using export GMX_NOOPTIMIZEDKERNELS=1 prior to running mdrun solves the issue for 4.5.4 but not 4.0.7 (on my architecture). To solve the problem in v4.0.7 (and this also works for v4.5.4) I must export GMX_NB_GENERIC=1; export GMX_NO_SOLV_OPT=1 prior to running mdrun. Note that the generic nonbonded kernel runs 1.1x slower then the unoptimized kernel, which itself runs 2.4x slower than the optimized kernel (on my architecture for this system). Thank you, Chris. -- original message -- Dear Users: I am trying to work with tabulated potentials for the first time. I would like to allow 2 charges to be at the same spot. Imagine taking a sodium and chloride ion and removing the LJ entirely and I want the system to evaluate energies correctly and be stable during the simulation. It was my intention to add a sort-of soft-core to the charge-charge interactions at very close distances so that the system did not crash. However, my testing didn't even get that far. I find that the Coulomb(SR) energy = nan even when using tables that suggest it should be zero. I suppose that there is some other place in the code where this problem develops? While looking at share/gromacs/top/table6-12.xvg it seems that this should not lead to energies of nan, but it does. $ head share/gromacs/top/table6-9.xvg # # Table AB1: ndisp=6 nrep=9 # 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 2.00e-03 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 4.00e-03 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 6.00e-03 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 8.00e-03 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 1.00e-02 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 1.20e-02 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 0.00e+00 and therefore, according to manual 4.5.4 page 151 equation 6.23, the [(qi*qj)/(4*pi*epsilon)]*f(rij) coulomb SR value should equal zero. I tried a few things, like modifying the values so that they are constant from 0 to 0.04 nm in the table.xvg file, but this did not help. here is the .top file that I used along with ffoplsaa to test out this idea: #include ffoplsaa.itp #include na.itp [ system ] ; name tables [ molecules ] ; name number MOL_A 1 MOL_B 1 and here is the .itp file that I used along with ffoplsaa to test out this idea: [ atomtypes ] aaa AA 8 15.999400.0 A0.0 0.0 [ moleculetype ] ; molname nrexcl MOL_A 1 [ atoms ] ; idat type res nr residu name at name cg nr charge 1 aaa 1 MLA AA 11.0 [ moleculetype ] ; molname nrexcl MOL_B 1 [ atoms ] ; idat type res nr residu name at name cg nr charge 1 aaa 1 MLB BB 1-1.0 And the .gro file: title 2 1MLA AA1 0.000 0.000 0.000 1MLB BB1 0.000 0.000 0.000 5.0 5.0 5.0 ### And my .mdp file is nsteps = 0 ns_type = grid rlist = 1 rcoulomb = 1 rvdw = 1 coulombtype = user tc_grps = System tau_t = 1.0 ld_seed = -1 ref_t = 300 gen_temp = 300 gen_vel = yes unconstrained_start = no gen_seed = -1 Pcoupl = berendsen pcoupltype = isotropic tau_p = 4 compressibility = 4.5e-5 ref_p = 1.0 Still, when I run a 0-step mdrun (mdrun -deffnm md1 -table mod2.table.xvg), I get: Energies (kJ/mol) LJ (SR) Coulomb (SR) PotentialKinetic En. Total Energy 0.0e+00nannannannan Temperature Pressure (bar) nannan # I get the same thing if the two charges are like or dislike... But if I then modity my topologies to turn the +1 charge to zero, I get sensible output:
[gmx-users] Rigid structure and flexible structure present in sing system
If I read between the lines correctly, you know how to do this in gromacs, but you wish that you got a big speedup from freezing most of the atoms in your system. If that is the case, then I think that gromacs can not help you in its current form. Therefore, I suggest that you try the Charmm software. It may be slower overall, but it gives you the speedup you expect when you freeze 95% of the atoms in your system so it may be much faster for your usage. Other software may also do this, but I have no experience with it. If you have explicit water, this won't be of much use since you're still only freezing 10% of the atoms in your system. If you simply want to know how to freeze one structure, then use freeze groups and energy exclusions without pressure coupling, or use position restraints and refcoordscaling=com with position restraints, or create an elastic network of restraints. Chris. -- original message -- Hi, For a while I've been looking for a way to include a rigid (macromolecular) structure and a flexible (small-protein) structure in a single system and have not had much success. I looked for a way to apply a different set of constraints to each structure, but couldn't find one. I looked into modifying the topology file, but the only method I could find was to artificially increase the mass of every atom in the rigid structure and this does not save any computing time. Does anyone know how I can fix the relative co-ordinates of one structure and only calculate the MD of another? Any constructive criticism (even if to inform me that I'm doing something ridiculous) is very gratefully received. Kind regards, Luke -- 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] beware gen-vel=no
Regarding gen-vel=no as discussed here: http://lists.gromacs.org/pipermail/gmx-users/2011-May/061151.html I would caution against the general use of gen-vel=no and then coupling to a temperature of 300 K with the Berendsen thermostat. One problem that can arise is concerted unfolding. Temperature is random velocity, but with gen-vel=no and a structure that is just slightly too compact, you will have a slight net force that gets massively scaled up by the coupling algorithm and results, not in true temperature, but in large scale coordinated motion. It's just something to look out for. gen-vel=no should not be a problem with the sd algorithm or with some type of position restrained equilibration procedure or, of course, if loading in velocities from some cpt or trr file. Chris. -- 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] Count mismatch for state entry SDx, code count is 754728, file count is 0
Dear Users: Using gromacs 4.0.5, I find that there are at least some cases where some type of disk error can get propagated through both my.tpr and my_prev.tpr, complicating restarts. This used to be a bigger problem in gromacs 3, and I don't recall ever seeing it in gromacs 4 so I thought I would post a notification. I'm just going to extract some coordinates and restart, but ideally this wouldn't happen. A google search for the relevant error Count mismatch for state entry only turns up some online source code. I don't know if this error occurs in 4.5.3, and it's not binary reproducible so that would be difficult to check. Still, the error checking that regularly occurs prior to overwriting the previous (and without error) _prev.cpt file with a new (and with error) _prev.cpt file seemed to not catch this problem, at least with gromacs 4.0.5. The run that wrote out the .tpr finished normally due to -maxh, with a stderr that looked like this: ... snip ... starting mdrun 'Generated by genbox' 1000 steps, 2.0 ps (continuing from step 3769350, 7538.7 ps). [gpc-f138n034:06165] 15 more processes have sent help message help-mpi-btl-base.txt / btl:no-nics [gpc-f138n034:06165] Set MCA parameter orte_base_help_aggregate to 0 to see all help / error messages Step 5036590: Run time exceeded 47.322 hours, will terminate the run Step 5036600: Run time exceeded 47.322 hours, will terminate the run Average load imbalance: 0.2 % Part of the total run time spent waiting due to load imbalance: 0.2 % Steps where the load balancing was limited by -rdd, -rcon and/or -dds: X 0 % Z 0 % Average PME mesh/force load: 0.745 Part of the total run time spent waiting due to PP/PME imbalance: 4.9 % Parallel run - timing based on wallclock. NODE (s) Real (s) (%) Time: 170485.000 170485.000100.0 1d23h21:25 (Mnbf/s) (GFlops) (ns/day) (hour/ns) Performance:625.583 31.889 1.284 18.685 gcq#165: I'm a Jerk (F. Black) gcq#165: I'm a Jerk (F. Black) # And then when I gmxcheck both of the .cpt files I get the exact same error, although the files do differ: $ diff md1.cpt md1_prev.cpt Binary files md1.cpt and md1_prev.cpt differ $ gmxcheck -f md1.cpt :-) G R O M A C S (-: S C A M O R G :-) 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. :-) gmxcheck (-: Option Filename Type Description -fmd1.cpt Input, Opt! Trajectory: xtc trr trj gro g96 pdb cpt -f2 traj.xtc Input, Opt. Trajectory: xtc trr trj gro g96 pdb cpt -s1 top1.tpr Input, Opt. Run input file: tpr tpb tpa -s2 top2.tpr Input, Opt. Run input file: tpr tpb tpa -c topol.tpr Input, Opt. Structure+mass(db): tpr tpb tpa gro g96 pdb -e ener.edr Input, Opt. Energy file: edr ene -e2 ener2.edr Input, Opt. Energy file: edr ene -n index.ndx Input, Opt. Index file -mdoc.tex Output, Opt. LaTeX file Option Type Value Description -- -[no]h bool no Print help info and quit -niceint0 Set the nicelevel -vdwfac real 0.8 Fraction of sum of VdW radii used as warning cutoff -bonlo real 0.4 Min. fract. of sum of VdW radii for bonded atoms -bonhi real 0.7 Max. fract. of sum of VdW radii for bonded atoms -tol real 0.001 Relative tolerance for comparing real values defined as 2*(a-b)/(|a|+|b|) -[no]ab bool no Compare the A and B topology from one file -lastenerstring Last energy term to compare (if not given all are tested). It makes sense to go up until the Pressure. Checking file md1.cpt --- Program gmxcheck, VERSION 4.0.5 Source code file: checkpoint.c, line: 186 Fatal error: Count mismatch for state entry SDx, code count is 754728, file count is 0 --- Confirmed (Star Trek) and the
[gmx-users] Re: Count mismatch for state entry SDx, code count is 754728, file count is 0
Apologies: my.tpr and my_prev.tpr should have read my.cpt and my_prev.cpt. On 11-05-05 12:36 PM, Chris Neale wrote: Dear Users: Using gromacs 4.0.5, I find that there are at least some cases where some type of disk error can get propagated through both my.tpr and my_prev.tpr, complicating restarts. This used to be a bigger problem in gromacs 3, and I don't recall ever seeing it in gromacs 4 so I thought I would post a notification. I'm just going to extract some coordinates and restart, but ideally this wouldn't happen. A google search for the relevant error Count mismatch for state entry only turns up some online source code. I don't know if this error occurs in 4.5.3, and it's not binary reproducible so that would be difficult to check. Still, the error checking that regularly occurs prior to overwriting the previous (and without error) _prev.cpt file with a new (and with error) _prev.cpt file seemed to not catch this problem, at least with gromacs 4.0.5. The run that wrote out the .tpr finished normally due to -maxh, with a stderr that looked like this: ... snip ... starting mdrun 'Generated by genbox' 1000 steps, 2.0 ps (continuing from step 3769350, 7538.7 ps). [gpc-f138n034:06165] 15 more processes have sent help message help-mpi-btl-base.txt / btl:no-nics [gpc-f138n034:06165] Set MCA parameter orte_base_help_aggregate to 0 to see all help / error messages Step 5036590: Run time exceeded 47.322 hours, will terminate the run Step 5036600: Run time exceeded 47.322 hours, will terminate the run Average load imbalance: 0.2 % Part of the total run time spent waiting due to load imbalance: 0.2 % Steps where the load balancing was limited by -rdd, -rcon and/or -dds: X 0 % Z 0 % Average PME mesh/force load: 0.745 Part of the total run time spent waiting due to PP/PME imbalance: 4.9 % Parallel run - timing based on wallclock. NODE (s) Real (s) (%) Time: 170485.000 170485.000100.0 1d23h21:25 (Mnbf/s) (GFlops) (ns/day) (hour/ns) Performance:625.583 31.889 1.284 18.685 gcq#165: I'm a Jerk (F. Black) gcq#165: I'm a Jerk (F. Black) # And then when I gmxcheck both of the .cpt files I get the exact same error, although the files do differ: $ diff md1.cpt md1_prev.cpt Binary files md1.cpt and md1_prev.cpt differ $ gmxcheck -f md1.cpt :-) G R O M A C S (-: S C A M O R G :-) 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. :-) gmxcheck (-: Option Filename Type Description -fmd1.cpt Input, Opt! Trajectory: xtc trr trj gro g96 pdb cpt -f2 traj.xtc Input, Opt. Trajectory: xtc trr trj gro g96 pdb cpt -s1 top1.tpr Input, Opt. Run input file: tpr tpb tpa -s2 top2.tpr Input, Opt. Run input file: tpr tpb tpa -c topol.tpr Input, Opt. Structure+mass(db): tpr tpb tpa gro g96 pdb -e ener.edr Input, Opt. Energy file: edr ene -e2 ener2.edr Input, Opt. Energy file: edr ene -n index.ndx Input, Opt. Index file -mdoc.tex Output, Opt. LaTeX file Option Type Value Description -- -[no]h bool no Print help info and quit -niceint0 Set the nicelevel -vdwfac real 0.8 Fraction of sum of VdW radii used as warning cutoff -bonlo real 0.4 Min. fract. of sum of VdW radii for bonded atoms -bonhi real 0.7 Max. fract. of sum of VdW radii for bonded atoms -tol real 0.001 Relative tolerance for comparing real values defined as 2*(a-b)/(|a|+|b|) -[no]ab bool no Compare the A and B topology from one file -lastenerstring Last energy term to compare (if not given all are tested). It makes sense to go up until the Pressure. Checking file md1.cpt --- Program gmxcheck, VERSION 4.0.5 Source code file: checkpoint.c, line: 186 Fatal error: Count mismatch for state entry SDx, code count is 754728, file count
[gmx-users] constraints
generally, look at mdout.mdp from grompp -- original message -- In my input file if I don't specify constraints then What is default constraints=none constraints=all bonds NIlesh -- 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] free energy
Your windows are restrained. The PMF that you get out of WHAM is a representation of the relative sampling after removing the umbrella biases. Sounds like you are saying that you look at the still-biased trajectories and you see different a different distribution of states than you do in the de-biased PMF. not sure what the problem is here. Perhaps go back to some review literature on US. Now, if you saw more bound than unbound in unrestrained simulations, then that's a different story, but that doesn't appear to be what you are talking about. Chris. -- original message -- Hi all I have performed a PMF simulation of taking part of a molecule out of the cavity of a host using umbrella sampling. The free energy curve suggests that the guest prefers to be outside the host as the bound state is higher in energy, or the free energy difference to go in is positive. However when I view the trajectories for each window it appears that there is always more bound states than unbound. This leads me to doubt my free energy profile? Cheers Gavin -- 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] free energy
yes it will, and so will it affect the profile if water molecules or ions go in when both chains are absent. You'll need to determine what question you are trying to answer and also think pretty hard about what your PMF really means in the context of this system. Chris. -- original message -- Hi Chris My windows are restrained obviously using the force constant in the mdp file. The trajectories that I have viewed are those of the individual biased sampling windows. I have not put on the unbiased simulations yet. There is also the following issue: The simulations involve two identical molecules containing hydrocarbon chains. I calculate the PMF to take a specific hydrocarbon chain of one molecule out of a specific site on the neighbouring molecule. When this guest chain goes out, other chains can go in (intramolecular or intermolecular). Will this affect the profile? Gavin -- 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] free energy
Dear Gavin: I'm not seeing a lot of effort on your part to answer these questions. I suspect that you can answer some of this. Good luck! Chris. -- original message -- The current simulations are currently in vacuum. Does the following mdp file seem ok. Note that this is a production run using the final configuration after a lengthy equilibration. Also I was thinking about trying to prevent the other chains from entering the specific site; is there a a way to implement this in gromacs. Gavin chris.neale at utoronto.ca wrote: yes it will, and so will it affect the profile if water molecules or ions go in when both chains are absent. You'll need to determine what question you are trying to answer and also think pretty hard about what your PMF really means in the context of this system. Chris. -- original message -- Hi Chris My windows are restrained obviously using the force constant in the mdp file. The trajectories that I have viewed are those of the individual biased sampling windows. I have not put on the unbiased simulations yet. There is also the following issue: The simulations involve two identical molecules containing hydrocarbon chains. I calculate the PMF to take a specific hydrocarbon chain of one molecule out of a specific site on the neighbouring molecule. When this guest chain goes out, other chains can go in (intramolecular or intermolecular). Will this affect the profile? Gavin -- 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
