[gmx-users] gmx 4.6.1, Expanded ensemble: weird balancing factors
Dear Gromacs users, I am trying to compute the free energy of hydration of ethane in order to get more familiar with the implementation of expanded ensemble in Gromacs v. 4.6.1. The simulations runs fine and the balancing/weight factors are equilibrated after around 1.6 ns (a short equilibration time but this should be ok for a rough estimate of dG). After this the simulation is continued for another 14 ns. At the end I look at the balancing factors to get a fast estimate of dG and see that they are tiny: MC-lambda information N CoulL VdwLCount G(in kT) dG(in kT) 1 0.000 0.000162900.0 -0.09268 2 0.200 0.00014940 -0.09268 -0.01647 3 0.500 0.00014627 -0.10914 -0.03867 << 4 1.000 0.00014080 -0.147810.02986 5 1.000 0.20014309 -0.117950.07408 6 1.000 0.40015170 -0.043860.01711 7 1.000 0.60015276 -0.02675 -0.02304 8 1.000 0.80015061 -0.04979 -0.06110 9 1.000 1.00014328 -0.110900.0 When adding these to get dG they value is way off what I get if I use the more traditional approach of separate simulations at each lambda value and using BAR at the end. Does anyone have any ideas why I get these values? Here is my mdp-file: integrator = sd tinit= 0 dt = 0.001 nsteps = 1500 comm-mode = Linear nstcomm = 10 nstlog = 1000 nstenergy= 100 nstlist= 10 ns_type= grid pbc= xyz rlist = 1.0 cutoff-scheme = verlet coulombtype= PME rcoulomb = 1.0 vdw-type = cutoff rvdw = 1.0 fourierspacing = 0.12 pme_order= 4 ewald_rtol = 1e-04 DispCorr= EnerPres tc-grps = System tau_t= 5.0 ref_t = 300 Pcoupl= berendsen tau_p = 10.0 compressibility = 4.5e-5 ref_p= 1.013 gen_vel = no constraints = all-bonds constraint-algorithm = lincs ; Free energy/expanded ensemble free-energy= expanded sc-alpha = 0.5 sc-r-power = 6 init-lambda-state = 0 coul-lambdas = 0.0 0.2 0.5 1.0 1.0 1.0 1.0 1.0 1.0 vdw-lambdas = 0.0 0.0 0.0 0.0 0.2 0.4 0.6 0.8 1.0 calc-lambda-neighbors = 9 nstdhdl = 200 dhdl-print-energy = yes couple-moltype = C1X couple-lambda0 = vdw-q couple-lambda1 = none couple-intramol= no nstexpanded = 100 lmc-stats= wang-landau lmc-move = metropolis lmc-weights-equil = wl-delta weight-equil-wl-delta = 0.001 wl-scale = 0.7 wl-ratio = 0.8 init-wl-delta = 1.0 wl-oneovert = yes Thanks in advance! / Joakim-- 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: Hamiltonian replica exchange umbrella sampling with gmx 4.6
Dear Michael, Thank you for your reply. Yes, it is relatively clear now. Will try to play around with this later today. Best, Joakim > Date: Thu, 7 Mar 2013 08:55:31 -0500 > From: Michael Shirts > Subject: Re: [gmx-users] Hamiltonian replica exchange umbrella > sampling with gmx 4.6 > To: Discussion list for GROMACS users > > Hi, Joakim- > > Hamiltonian exchange only should work if there is a lambda coupling > parameter that defines the potential at each state. You need to > define your pulling potential so that the coupling-lambda parameter > can be used to define the different pulling location centers along > your trajectory. Does this make it clearer? -- 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] Hamiltonian replica exchange umbrella sampling with gmx 4.6
Dear gmx-users, I am currently trying to runt Hamiltonian replica exchange umbrella sampling in hope to do some better sampling. I have generated a number of tpr-files along my reaction coordinate and they all run fine in independent simulations. The issue comes when I would like the replica exchange to start. The following line is used to initiate the exchange: mdrun_mpi -replex 1000 -s md -pf pullf -px pullx -multi 40 -maxh 0.5 All replicas have the same temperature and the following error is what I face seconds after submitting the job: The properties of the 40 systems are all the same, there is nothing to exchange For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors I could simply change the temperature between the replicas by 0.001 K and it would run I think. But that is not very elegant. Does anyone have any suggestions? Thanks in advance! Best regards, Joakim-- 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: The sum of the two largest charge group radii is larger than rlist - rvdw (because rlist < rvdw)
Dear Chris, We have used the following settings in our recent work (http://pubs.rsc.org/en/content/articlelanding/2013/CP/C3CP44472D): coulombtype=pme rcoulomb=1.0 rlist=1.0 vdw-type=cut-off rvdw=1.4 dispcorr=enerpres And this works fine without any notes etc from grompp, and you do not have to use the switching function. Should probably change the md.mdp file on the website...thanks for bringing this to our attention! A colleague of mine have used Slipids with a shorter rvdw and the bilayer properties were still very stable. When I use rcoulomb=rvdw=1.0 I see no major difference either except for a relatively big boost in performance. But please double check that your simulations are able to reproduce experimental observables if changing rvdw to around 1.0 nm., but it should be fine. Best regards, Joakim Date: Tue, 26 Feb 2013 19:56:12 + From: Christopher Neale Subject: [gmx-users] The sum of the two largest charge group radii is larger than rlist - rvdw (because rlist < rvdw) To: "gmx-users@gromacs.org" > Dear users: > > I am experimenting with the "Stockholm" lipid parameters (Slipids). I > downloaded the recommended .mdp file from the developers of this force field > (http://people.su.se/~jjm/Stockholm_Lipids/Downloads_files/md.mdp ) and was > surprised to see the treatment of non-bonded interactions: > > coulombtype = pme > rcoulomb = 1.0 > rlist = 1.0 > rlistlong = 1.6 > rvdw = 1.5 > rvdw-switch = 1.4 > vdw-type = switch > > it struck me as strange to have rlist < rvdw (I didn't even know that this > was allowed, but I tried it and it is) > > (Note that this is from their website, but their original paper would have > used rcoulomb = 1.4 and rlist = 1.4). > > When I run grompp version 4.5.5 with these parameters, I get the following > note: > > NOTE 2 [file empty.mdp]: > The sum of the two largest charge group radii (0.00) is larger than > rlist (1.00) - rvdw (1.50) > > This system has only one Na and one Cl ion (chosen for simplicity during mdp > testing). By moving these ions apart and doing a 0-step mdrun, I was able to > verify that the LJ interaction energy moves from the LJ(SR) to the LJ(LR) > component for d>=1 nm and that this energy is non-zero until d>=1.5 nm. > > I am asking if anybody sees a problem because (a) grompp threw a note, and > (b) even if this works alright for standard simulations, I just want to be > sure that if I get into more exotic parts of the gromacs code then I am not > likely to encounter situations in which this cutoff scheme introduces > problems (silent or otherwise). > > 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] Expanded Ensemble and Gromacs 4.6
Hi, I am trying to use the Expanded Ensemble (EE) method to compute the free energy of solvation of a small organic molecule. Basically I am playing around but I cannot get the simulations to run. Here are my EE and free energy settings: %--- free-energy = expanded ; Expanded ensemble couple-moltype = C1X ; Molecule to introduce couple-lambda0 = none ; Go from no interactions with solvent... couple-lambda1 = vdw-q ; ...to full interactions couple-intramol = no; Do not decouple internal interactions init_lambda_state = 0 ; Start from the first column of the lambda vectors delta-lambda= 0 ; No increments in lambda nstdhdl = 100 ; Frequency for writing dH/dlambda coul-lambdas= 0.0 0.25 0.5 0.75 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0 vdw-lambdas = 0.0 0.0 0.0 0.0 0.0 0.2 0.4 0.6 0.7 0.8 0.9 1.0 dhdl-print-energy = yes ; Include total energy in the dhdl file sc-alpha= 0.5 ; Soft core alpha parameter sc-power= 1 ; Power of lambda in soft core function sc-r-power = 6 ; Power of radial term in the soft core function sc-sigma= 0.3 ; Soft core sigma sc-coul = no; No soft core for Coulomb separate-dhdl-file = yes ; Seperate dhdl files dhdl-derivatives= yes ; Print derivatives of the Hamiltonian ; expanded ensemble variables nstexpanded = 100 ; Number of steps between attempts to change the Hamiltonian lmc-stats = wang-landau ; WL algorithm to explore state space lmc-move= gibbs ; Decides which state to move to lmc-weights-equil = number-steps ; EE weight updating stops after... weight-equil-number-steps = 500 ; ...10ns of simulation ; Seed for Monte Carlo in lambda space lmc-seed= 1993 lmc-repeats = 1 lmc-gibbsdelta = -1 lmc-forced-nstart = 0 symmetrized-transition-matrix = no nst-transition-matrix = -1 mininum-var-min = 100 wl-scale= 0.6 wl-ratio= 0.8 init-wl-delta = 1 wl-oneovert = yes %--- Besides this I use LINCS to constrain all bonds, sd-integrator and a "normal" cut-off scheme with a 1 fs time step. However, once I try to run the files on the cluster I always end up with LINCS warnings and after a few seconds the program crashes due to too many LINCS warnings. If I increase the time step to 2 fs I run into problems SETTLE for the water. I always start from a nice structure that is taken as the final snap shot from a simple 1 ns MD simulation of my system. What could be the problem? If I use free_energy = yes instead of EE things work fine. Did I perhaps mess up the EE settings or something? Thanks in advance. Best, Joakim -- 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
Re: [gmx-users] Pre-equilibrated CHARMM lipid bilayers
We have performed relatively long simulations of a WALP23 peptide embedded in DLPC and DOPC bilayers with different flavors of the Amber FF family. Currently I am working with much bigger protein and see good agreement between the simulations and experiments. You could probably use those bilayers with CHARMM, just make sure you equilibrate a while before just to be safe. Best, Joakim On Dec 19, 2012, at 19:22 , gmx-users-requ...@gromacs.org wrote: > Date: Wed, 19 Dec 2012 18:45:33 +0300 > From: James Starlight > Subject: Re: [gmx-users] Pre-equilibrated CHARMM lipid bilayers > To: Discussion list for GROMACS users > Message-ID: > > Content-Type: text/plain; charset=ISO-8859-1 > > Hi Joakim! > > MAny thanks for such big collection of the bilayers. Have your > simulated different protein embedded in that bilayers ? What force > field did you use for parametrisation of the proteins for such > simulations ? From the mdp file I noticed relatively big cutoffs for > vdw as well as electrostatics ( bigger than in charmm force fields). > Would it possible to simulate proteins parametrized by charmm in that > lipids ussing charmm parametres for whole system? > > > Thanks again, > > 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
Re: [gmx-users] Pre-equilibrated CHARMM lipid bilayers
Hi James, You could also have look at our website: http://people.su.se/~jjm There you find a number of bilayers ready to use (among them there is a POPE bilayer). They were generated with another FF but if you wish to use CHARMM you can since the atom numbering/naming is the same for Slipids and CHARMM. Best, Joakim Message: 2 Date: Tue, 18 Dec 2012 21:07:22 -0800 From: James Starlight Subject: Re: [gmx-users] Pre-equilibrated CHARMM lipid bilayers To: Discussion list for GROMACS users Message-ID: Content-Type: text/plain; charset=ISO-8859-1 Justin, thanks again. As I understood gromacs already had had parameters for charmm lipid so the main approach is to do ITP file for 1 lipid by means of pdb2gmx isnt it? By the way is there any way to convert PSF or CRD file to PDB? I've found suitable bilayer for my simulation but it lack such coordinates. POPE Bilayer (310.15K, A=65.2 Â/lipid, 10ns, 340 lipids, noLRCs): CHARMM PSF, X-Plor/NAMD PSF, and CHARMM CRD 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] forcefields for lipids
Hi Tom, we have parameters for unsaturated tails now as well as for PE lipids ready for use, the manuscript for these have been submitted. PG and PS lipids are on the way and parameters for cholesterol and sphingomyelin are finished now and the manuscript should be finished relatively soon. Best, J On May 22, 2012, at 15:40 , gmx-users-requ...@gromacs.org wrote: > Date: Tue, 22 May 2012 13:49:48 +0100 > From: Thomas Piggot > Subject: Re: [gmx-users] forcefields for lipids > To: Discussion list for GROMACS users > Message-ID: <4fbb8b6c.70...@soton.ac.uk> > Content-Type: text/plain; charset="ISO-8859-1"; format=flowed > > Hi Joakim, > > I am interested to know what other lipids this force field has been > extended to include. Do you have parameters for unsaturated lipid tails > and other classes of phospholipid (such as PE and PG)? > > Cheers > > Tom -- 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] forcefields for lipids
Hi, have a look at out recently developed all-atom force field. It is compatible with AMBER force fields for proteins (ff99SB, ff99SB-ILDN and ff03). The compatibility was tested via microsecond simulations of a peptide embedded in a bilayer. The paper describing the peptide simulations along with an extension of the parameter set has been submitted, but for now have a look at this published paper: J. Phys. Chem. B, 2012, 116 (10), 3164-3179, doi: 10.1021/jp212503e All our parameters are available at on Lipidbook. Let me know if you need anything else! Best, J -- 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
