[gmx-users] Determining Force Constants for CG modelling
Hello community, I'm trying to create a topology for a molecule using the MARTINI force field which is a coarse-grain (CG) forcefield. I understand that to optimize the bonded parameters, one needs to model the AA version and extract the equilibrium angle and force constants. As is said in the MARTINI manual, an example is "the angle distribution function of a CG triplet can be directly compared to the distribution function obtained from the AA simulation of the centers of mass of the corresponding atoms". Once you find the pseudo-CG beads from the AA model, I see how you would get the equilibrium parameters, but not the force constants. Is it simply taking the bonded potentials from the AA simulation and have them matched to the CG bonded potentials by changing the force constants so they have the same potential? Is it necessary to use the Boltzmann Inversion Approach? Thank you for your help! -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering -- gmx-users mailing listgmx-users@gromacs.org http://lists.gromacs.org/mailman/listinfo/gmx-users * Only plain text messages are allowed! * 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: Chemical Potential
Thank you Jan. -Fabian On Tue, May 22, 2012 at 6:54 PM, Fabian Casteblanco wrote: > Hello community, > > I'm just trying to explore what kind of calculations one can do on > polymer systems (pure or in water) in order to validate the force > field works accurately for that system. I know there are basics such > as density, volume, dH of vaporization, isothermal compressibility, > heat capacity, etc. I've been reading about the particle insertion > method to calculate chemical potentials. Since the chemical potential > is simply the change in gibbs as the number of particles changes, can > one use the g_bar method to simply insert/delete a molecule to/from > the system? > > Anybody know an article where this or something similar was done? Thanks. > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering -- 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] Chemical Potential
Hello community, I'm just trying to explore what kind of calculations one can do on polymer systems (pure or in water) in order to validate the force field works accurately for that system. I know there are basics such as density, volume, dH of vaporization, isothermal compressibility, heat capacity, etc. I've been reading about the particle insertion method to calculate chemical potentials. Since the chemical potential is simply the change in gibbs as the number of particles changes, can one use the g_bar method to simply insert/delete a molecule to/from the system? Anybody know an article where this or something similar was done? Thanks. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering -- 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] Slow-growth method question
Hello community, I recently ran several simulations using the g_bar method on several molecule decouplings and I even applied it to some functional group mutations and it seems to work out quite well so far. The only thing is that it requires several simulations to run. I'm attempting to replicate this group mutation using the original slow-growth method that others have done in the past just to see how well it compares to the g_bar method. I can't seem to find any sort of instructions on running this method on Gromacs.I have the parameters from g_bar below but I'm trying to figure out how to alternate them below for the slow growth method. Is it simply the delta_lambda parameter that needs to be changed? Would I have to divide the timesteps by delta_lambda to make sure it goes up to Lambda=1? How do you calculate the end result?Thank you. ; Free energy control stuff free_energy = yes ;*** - Indicates we are doing a free energy calculation, and that interpolation between the A and B states of the chosen molecule (defined below) will occur. init_lambda = 0.0 ;*** - Value of λ delta_lambda= 0 ;*** - The value of λ can be incremented by some amount per timestep (i.e., δλ/δt) in a technique called "slow growth." This method can have significant errors associated with it, and thus we will make no time-dependent changes to our λ values. foreign_lambda = 0.05 ;*** - Additional values of λ for which ΔH will be written to dhdl.xvg (with frequency nstdhdl). The configurations generated in the trajectory at λ = init_lambda will have ΔH calculated for these same configurations at all values of λ = foreign_lambda. sc-alpha= 0.5 ;*** - the soft-core parameter, a value of 0 results in linear interpolation of the LJ and Coulomb interactions sc-power= 1.0 ;*** - the power for lambda in the soft-core function, only the values 1 and 2 are supported sc-sigma= 0.3 ;*** - the soft-core sigma for particles which have a C6 or C12 parameter equal to zero or a sigma smaller than sc_sigma ;couple-moltype = SIM ;*** - name of moleculetype to decouple ;couple-lambda0 = vdw ;*** - all interactions are on at lambda=0 ;couple-lambda1 = none ;*** - only vdw interactions are on at lambda=1 couple-intramol = no ;*** - Do not decouple intramolecular interactions. That is, the λ factor is applied to only solute-solvent nonbonded interactions and not solute-solute nonbonded interactions. nstdhdl = 10 ;*** - The frequency with which ∂H/∂λ and ΔH are written to dhdl.xvg. A value of 100 would probably suffice, since the resulting values will be highly correlated and the files will get very large. You may wish to increase this value to 100 for your own work. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering -- 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] FEP
Hello all, This is in reply to Michael shirts a while ago on a FEP of a R-CH3 to an R-H group. Below is the orignal email. I recently tested out mutating a CH3-CH3-(3dummy atoms) molecule on both sides in order to test out that a peturbation would give you a total of 0. forcefield used was CGenFF So for procedure was that I turned the CH3 group on the left to an H with 3 dummy atoms and I turned an H and 3 dummy atoms on the right to a CH3 group, essentially mutating the molecule to another version of itself (should give you a total dG of 0 if it works). STEP 1. (uncharging everything except the -CH2 group inbetween both sides CH3-CH2-HD3) *D are dummy atoms [ atoms ] ; nr type resnr residatom cgnr chargemasstotal_charge 1CG331 1 SIM C1 1-0.27 12.0110 CG331 0.00 12.0110 2 HGA3 1 SIM H11 2 0.09 1.00800 HGA30.00 1.00800 3 HGA3 1 SIM H12 3 0.09 1.00800 HGA30.00 1.00800 4 HGA3 1 SIM H13 4 0.09 1.00800 HGA30.00 1.00800 5CG331 1 SIM C2 5-0.27 12.0110 CG331 -0.1812.0110 6 HGA3 1 SIM H21 6 0.09 7 HGA3 1 SIM H22 7 0.09 8 HGA3 1 SIM H23 8 0.09 1.00800 HGA30.00 1.00800 9 DUM 1 SIM H31 9 0.00 1.00800 10 DUM 1 SIM H3210 0.00 1.00800 11 DUM 1 SIM H3311 0.00 1.00800 STEP 2. (mutation of groups. CH3-CH2-HD3 becomes D3H-CH2-CH3) [ atoms ] ; nr type resnr residatom cgnr chargemasstotal_charge 1CG331 1 SIM C1 1 0.00 12.0110 HGA30.00 1.00800 2 HGA3 1 SIM H11 2 0.00 1.00800 DUM 0.00 1.00800 3 HGA3 1 SIM H12 3 0.00 1.00800 DUM 0.00 1.00800 4 HGA3 1 SIM H13 4 0.00 1.00800 DUM 0.001.00800 5CG331 1 SIM C2 5 -0.18 6 HGA31 SIM H21 6 0.09 7 HGA31 SIM H22 7 0.09 8 HGA31 SIM H23 8 0.00 1.00800 CG331 0.00 12.0110 9 DUM1 SIM H31 9 0.00 1.00800 HGA30.00 1.00800 10 DUM1 SIM H3210 0.00 1.00800 HGA3 0.00 1.00800 11 DUM1 SIM H3311 0.00 1.00800 HGA30.00 1.00800 STEP3. (opposite of step1) [ atoms ] ; nr type resnr residatom cgnr chargemasstotal_charge 1 HGA3 1 SIM C1 1 0.00 12.0110 HGA3 0.091.00800 2 DUM 1 SIMH11 2 0.00 1.00800 DUM0.00 1.00800 3 DUM 1 SIMH12 3 0.00 1.00800 DUM0.00 1.00800 4 DUM 1 SIMH13 4 0.00 1.00800 DUM0.00 1.00800 5CG331 1 SIM C2 5-0.18 12.0110 CG331 -0.27 12.0110 6 HGA3 1 SIM H21 6 0.09 7 HGA3 1 SIM H22 7 0.09 8CG331 1 SIM H23 8 0.00 1.00800 CG331 -0.27 12.0110 9 HGA3 1 SIM H31 9 0.00 1.00800 HGA30.09 1.00800 10 HGA3 1 SIM H3210 0.00 1.00800 HGA30.09 1.00800 11 HGA3 1 SIM H3311 0.00 1.00800 HGA30.09 1.00800 So in the end, the procedure worked using the g_bar method. I took 20 simulations for steps 1 and 3 and 50 simulations for step 2. Here are the results: Step1: -37.62 +/- 0 Step2: -0.12 +/- 0.24 Step3: 37.68 +/- 0.13 TOTAL ~ -0.06 kJ/mol which is expected My question is that I did something similar for a larger molecule (65 atoms) and someone said before that step 1 and 3 should be very small. I'm getting values in the near -20 kJ/mol for step 1 and 3. I feel it should be straightforward since I'm simply uncharging a CH3 group and I don't get why I'm getting such a large number. If it is suppose to be very small, is it safe to assume
[gmx-users] couple-moltype question
Hello Gromacs community, I am trying to simply decharge a part of large molecule. I know from the tutorial we can use 'couple-moltype' along with 'couple-lambda0', etc, but in this case I simply change the topology to state A and then I have state B written with no charges since I'm only doing a piece of the molecule. I have a free energy but it is only the free energy including interactions between the piece of the molecule and surrounding solvent. I see 'couple-intramol' is an option to have intramolecular interactions also decoupled but it requires a 'couple-moltype' to be set. Is there any way I can set that variable to simply a piece of a molecule? Thank you in advance for your help. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering -- 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] Radial distribution function by COM
Hello everyone, Is there any way to do a g(r) plot between the COM of a single solute particle, and the COMs of each solvent molecule around it? It seems to only let your choice be the COM for the first pick. Is there any way to do it for both choices? Thanks. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering -- 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: Determining energies between a solute and solvent
Thanks. This helps. One more question, is it possible to further break down the energy interactions between a piece of a protein and another piece of the solvent? For example, if I wanted the Cou and vdw interactions specifically for an -OH group on a solvent and some chain on the solute. I know how to do this using make_ndx for g(r) plots but I don't see how you can apply it to energygrps on the *.mdp file. Thanks again for all your help. -Fabian On Thu, Feb 23, 2012 at 2:41 PM, Fabian Casteblanco wrote: > Thanks. So before I run the simulation, I must use energygrps (for > example: energygrps: Protein SOL) in my md.mdp file so that it knows > to write down specific interactions between the two groups. Will > then it simply appear when I later run g_energy? > > -Fabian > > > - > Short-range nonbonded interactions can be decomposed using proper energygrps > in > the .mdp file. > > Perhaps, but g_enemat also requires that energygrps are specified when running > the simulation. Otherwise, the applicable terms in the .edr file are not > decomposed. > > -Justin > > > > On Thu, Feb 23, 2012 at 1:56 PM, Fabian Casteblanco > wrote: >> Hello everyone, >> >> Is it possible to see the energy (LJ, Cou) for simply the solute >> interacting with the solvent? >> >> For example, g_energy will calculate all the energies for the entire >> system interactions which include solvent-solvent interactions. I >> would simply like the solute-solvent interactions or possibly simply >> the solvent-solvent interaction energies. I see g_enemat as a >> possible function to use. Is this the best way? >> >> Thanks. >> >> -- >> Best regards, >> >> Fabian F. Casteblanco >> Rutgers University -- >> Chemical Engineering > > > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering -- 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: Determining energies between a solute and solvent
Thanks. So before I run the simulation, I must use energygrps (for example: energygrps: Protein SOL) in my md.mdp file so that it knows to write down specific interactions between the two groups. Will then it simply appear when I later run g_energy? -Fabian - Short-range nonbonded interactions can be decomposed using proper energygrps in the .mdp file. Perhaps, but g_enemat also requires that energygrps are specified when running the simulation. Otherwise, the applicable terms in the .edr file are not decomposed. -Justin On Thu, Feb 23, 2012 at 1:56 PM, Fabian Casteblanco wrote: > Hello everyone, > > Is it possible to see the energy (LJ, Cou) for simply the solute > interacting with the solvent? > > For example, g_energy will calculate all the energies for the entire > system interactions which include solvent-solvent interactions. I > would simply like the solute-solvent interactions or possibly simply > the solvent-solvent interaction energies. I see g_enemat as a > possible function to use. Is this the best way? > > Thanks. > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering -- 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] Determining energies between a solute and solvent
Hello everyone, Is it possible to see the energy (LJ, Cou) for simply the solute interacting with the solvent? For example, g_energy will calculate all the energies for the entire system interactions which include solvent-solvent interactions. I would simply like the solute-solvent interactions or possibly simply the solvent-solvent interaction energies. I see g_enemat as a possible function to use. Is this the best way? Thanks. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering -- 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] Noise in a radial distribution function using g_rdf
Sorry, I meant on the previous post, using g_rdf On Mon, Feb 13, 2012 at 12:23 PM, Fabian Casteblanco wrote: > Hello everyone, > > I'm using radial distribution functions on gromacs for the first time. > Is it normal to see some noise on these plots? You can still see a > trend on the plot but it appears somewhat noisy and not as nice and > linear as some other plots I've seen in other simulations like the > tutorial. > > Thanks. > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering -- 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] Noise in a radial distribution function using g_bar
Hello everyone, I'm using radial distribution functions on gromacs for the first time. Is it normal to see some noise on these plots? You can still see a trend on the plot but it appears somewhat noisy and not as nice and linear as some other plots I've seen in other simulations like the tutorial. Thanks. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering -- 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: Free Energy tutorial - choosing number of solvent molecules
Hello Justin, I am running 2,000,000 time steps for the actual MD run (4,000 ps). Depending on how many solvent molecules I start with, I get slightly different results. Do you think its ok to run several different tests and take the average? Or perhaps take the end results of a shorter MD run and use those as the starting coordinates for a new run? Thanks, Fabian Fabian Casteblanco wrote: > Hello all, > > I'm running the same process from the free energy tutorial by Justin > Lemkul...http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/index.html. > > How did the number of solvent particles get chosen (in the tutorial, 210 > molecules were chosen)? I seem to be getting slightly different If memory serves, I reproduced what was in the original paper, but do check. > results (ranging from as small as 200 J/mol to about 5 kJ/mol depending > on how many molecules I choose (ranging for example from 210 ethanol > molecules to about 610 ethanol molecules for the largest energy > difference change of about 5 kJ). I keep running tests to see if there > is some sort of minimum atom number to get steady consistent numbers but > I can't seem to find it. When I plot for example the bar.xvg & > barint.xvg for both sets to see where the lines don't match up, its > usually one or two points that differ slightly which cause the free > energies in the end to be slightly different.I seem to be noticing > too that the more atoms I use, the free energy gets a little bit lower. > > Does anybody have any experience with this? > How long are your simulations? I have experienced the case (using a water-octanol solvent mixture) where the initial configuration made a big difference in the result, so longer simulations and multiple configurations for the solvent were necessary to get reliable averages. -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 ======== On Thu, Jan 26, 2012 at 11:21 AM, Fabian Casteblanco wrote: > > Hello all, > > I'm running the same process from the free energy tutorial by Justin > Lemkul...http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/index.html. > > How did the number of solvent particles get chosen (in the tutorial, 210 > molecules were chosen)? I seem to be getting slightly different results > (ranging from as small as 200 J/mol to about 5 kJ/mol depending on how many > molecules I choose (ranging for example from 210 ethanol molecules to about > 610 ethanol molecules for the largest energy difference change of about 5 > kJ). I keep running tests to see if there is some sort of minimum atom > number to get steady consistent numbers but I can't seem to find it. When I > plot for example the bar.xvg & barint.xvg for both sets to see where the > lines don't match up, its usually one or two points that differ slightly > which cause the free energies in the end to be slightly different. I seem > to be noticing too that the more atoms I use, the free energy gets a little > bit lower. > > Does anybody have any experience with this? > > Thanks. > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > C: +908 917 0723 > E: fabian.castebla...@gmail.com > -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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 tutorial - choosing number of solvent molecules
Hello all, I'm running the same process from the free energy tutorial by Justin Lemkul... http://www.bevanlab.biochem.vt.edu/Pages/Personal/justin/gmx-tutorials/free_energy/index.html . How did the number of solvent particles get chosen (in the tutorial, 210 molecules were chosen)? I seem to be getting slightly different results (ranging from as small as 200 J/mol to about 5 kJ/mol depending on how many molecules I choose (ranging for example from 210 ethanol molecules to about 610 ethanol molecules for the largest energy difference change of about 5 kJ). I keep running tests to see if there is some sort of minimum atom number to get steady consistent numbers but I can't seem to find it. When I plot for example the bar.xvg & barint.xvg for both sets to see where the lines don't match up, its usually one or two points that differ slightly which cause the free energies in the end to be slightly different.I seem to be noticing too that the more atoms I use, the free energy gets a little bit lower. Does anybody have any experience with this? Thanks. -- *Best regards,* ** *Fabian F. Casteblanco* *Rutgers University -- * *C: +908 917 0723* *E: **fabian.castebla...@gmail.com* -- 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] decoupling a group of a molecule
Hello, Does anybody know how to do the same decoupling technique from Justin's tutorial but only for a small piece of a molecule rather than the entire molecule? Is it simply writing on the topology B state as zero charge and then mutating to a dummy atom? Thanks. -- Best regards, Fabian F. Casteblanco Rutgers University -- -- 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 of a mutated molecule
Hello all,
Please if anybody can help.
I'm trying to mutate a -CH3 to an -H (I guess with 3 dummy atoms
attached to it). Below I sketched the process. I broke it up into
3 steps and I wanted to use g_bar for the actual mutation step (step
2). Does anybody have any experience with something like this? I
keep getting LINCS errors like this one:
Step 3, time 0.006 (ps) LINCS WARNING
relative constraint deviation after LINCS:
rms 0.000378, max 0.010709 (between atoms 9 and 68)
bonds that rotated more than 30 degrees:
atom 1 atom 2 angle previous, current, constraint length
9 68 57.10. 0.1099 0.
9 66 90.00. 0.1122 0.
... for some of the Lambdas in the mutation part. I understand this
is a common error and it means that things are changing so my thinking
is that it is trying to rotate the dummy bonds for some reason I don't
understand. Are my steps below correct? Are dummy atoms suppose to
keep their original mass (not zero) ? Am I suppose to add parameters
for State B for these dummy atom bonds/angles also?
If anybody can help, I would greatly appreciate it. Thank you.
.--{Atom 68 (H --> Dum)}
|
R - {Atom 8 C} - {Atom 9 (C --> H)} - {Atom 66 (H --> Dum)}
|
'--{Atom 67 (H --> Dum)}
Essentially.. a R-C-[CH3] --> R-C-[H]
--
Force Field: CGenFF
STEP 1: Decoupling: -[CH3] (0 Charge, No LJ Interactions)
[ atoms ]
; nr type resnr resid atom cgnr charge
mass typeB chargeB massB
8 CG301 1 SIM C38 0.00
9 CG331 1 SIM C25 9 -0.27
12.01100 CG331 0.00 12.01100 ;
66 HGA3 1 SIM H251 66 0.09
1.00800 HGA3 0.00 1.00800 ;
67 HGA3 1 SIM H252 67 0.09
1.00800 HGA3 0.00 1.00800 ;
68 HGA3 1 SIM H253 68 0.09
1.00800 HGA3 0.00 1.00800 ;
--
STEP 2: Mutation: -[CH3] to -[H]
[ atoms ]
; nr type resnr residatomcgnrcharge
masstypeB chargeBmassB
8 CG3011SIM C38 0.00
12.01100
CG3110.0012.01100 ;
9 CG3311SIM C25 9 0.00 12.01100
HGA1 0.001.00800;
66 HGA3 1SIM H251660.00
1.00800
DUM 0.00 0.00 ;
67 HGA3 1SIM H252670.00 1.00800
DUM 0.00 0.00 ;
68 HGA3 1SIM H253680.00
1.00800DUM 0.00 0.00 ;
--
STEP 3: Coupling: -[H] (LJ Interactions, Full Charge)
[ atoms ]
; nr type resnr residatomcgnrcharge
masstypeB chargeBmassB
8 CG3111SIM C38 0.00 12.01100
CG311 -0.09 12.01100 ;
9 HGA1 1SIM C25 9 0.00
1.00800 HGA10.09 1.00800;
66 DUM 1SIM H25166 0.00 0.00
DUM0.00 0.00 ;
67 DUM 1SIM H25267 0.00 0.00
DUM 0.00 0.00 ;
68 DUM 1SIM H25368 0.00 0.00
DUM0.00 0.00 ;
--
Best regards,
Fabian F. Casteblanco
Rutgers University --
Chemical Engineering PhD Student
E: fabian.castebla...@gmail.com
--
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] Mixing Force Fields
Hello community, I have a general question about mixing force fields. I want to try to use CGenFF to model the drug and I wanted to use CHARMM for the solvent. I see the general advice given online is to always use only one forcefield but I know CGenFF was developed off CHARMM but specifically for drugs. If anything, I would like to try it to see how different results are from just using CGenFF. How would one go about using two forcefields on your topology? For example below, would one simply add the non-default FF to the default FF? ; File 'topol.top' was generated ; By user: Administrator (500) ; On host: fabian-3bef2d6a ; At date: Fri April 1st 02:59 PM 2011 ; ; This is your topology file ; DRUG SOLUTION w/ Ethanol ; Include forcefield parameters #include "Charmm.ff" #include "Cgenff.ff" **Thanks for everyones help. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Peturbing a Dihedral for FEP
Hello all, It seems I'm still getting errors when doing a FEP on a molecule (a -CH3 to a -H). This below is for when I was charging things from a -H uncharged to -H charged, although it also happens when I'm actually converting the -CH3 to -H (at Lambdas greater than 85%). I made sure to keep the charges balanced at 0 while mutating and I did it at 3 steps like Michael Shirts suggested. Set 1: turn R-CH3 charges off in a way that preserves the total charge. Set 2: change CH3 LJ to H LJ Set 3: Turn R-H charges on in a way that preserves the total charge. In the portion of the error below atom 9 is -C9-(H67,H68,H66) which in this specific case H67 is already a dummy molecule with no mass or charge. From what I can see, it seems that the atoms do not know how to treat the dummy molecules in terms of angles. How should I treat the dummy molecules? Should I be treating them like hollow spheres with no charge so I would assign them angle constraints? I think it can also be that I'm peturbing the dihedral angles incorrectly. I received errors at first saying that dihedral multiplicities can't be peturbed so I had to equal the multiplicities just to get it to run. Does anybody have any experience with this? Thank you for your help. -Fabian Casteblanco Portion of Error Output: - Reading file nvt0.5.tpr, VERSION 4.5.3 (single precision) starting mdrun 'SIMVASTATIN' 15 steps, 300.0 ps. Step 7, time 0.014 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.006111, max 0.139443 (between atoms 9 and 67) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length Step 8, time 0.016 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.007341, max 0.167622 (between atoms 9 and 67) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length Step 8, time 0.016 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.008771, max 0.201182 (between atoms 9 and 67) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length > On Mon, Oct 10, 2011 at 4:10 PM, Fabian Casteblanco > wrote: >> Hi all, >> >> I have an additional question related to what Steven Neumann was >> mentioning. I actually have to do a molecule mutation. I'm trying to >> use Michael Shirts method 1) making small >> changes 'alchemical' changes in the molecules and computing the free >> energies by any method (BAR, TI, etc). I'm specifically want to try >> to use BAR at the end once I collect all the data. This helped a lot >> on clarification since it seemed that Justin's tutorial is essentially >> a FEP except its using the BAR mathematical method for computing the >> complete decoupling of the molecule rather than using the old FEP >> mathematics of the exponential averaging formula. So BAR is only >> referring to the mathematical code used to calculate the overall free >> energy for the FEP, correct? >> >> My question is, for a mutation of a -CH3 group to a -H group, is it >> better to simply run: >> [+ from (Lambda=0 , R-CH3, full charges and interactions -STATE A) >> --> (Lambda=1, R-CH, full charges and interactions -STATE B)] >> >> OR >> >> [1) from (Lambda=0 , R-CH3, STATE A : Charges and LJ Interactions: ON) >> 2) (Lambda=?, -CH3 Charges: OFF ,LJ Interactions: ON and unmutated) >> 3) (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ Interactions: OFF) >> 4) (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ interactions: ON and Mutated to >> -H) >> 5) (Lambda=1, R-CH3, -CH3 STATE B : Charges and LJ Interactions: ON) >> >> Reason I'm asking is because when I try the first choice to do it >> STATE A to STATE B in one step, when I reach Lambda=0.85 and above on >> the NVT equilibration right after EM, I receive errors saying that >> bonds are moving way to far off their constraints which leads me to >> believe that the system is moving too far from where it was energy >> minimized. Errors such as: >> >> Step 188, time 0.376 (ps) >> LINCS WARNING >> relative constraint deviation after LINCS: >> rms 0.17, max 0.000636 (between atoms 9 and 68) >> bonds that rotated more than 30 degrees: >> atom 1 atom 2 angle previous, current, constraint length >> 9 68 31.2 0. 0.1110 0. >> >> Step 188, time 0.376 (ps) LINCS WARNING >> relative constraint deviation after LINCS: >> rms 0.15, max 0.000531 (between atoms 9 and 68) >> bonds that rotated more than 30 degrees: >> atom 1
[gmx-users] Re: FEP
Hello Michael or others, It seems I'm still getting errors when doing a FEP on a molecule (a -CH3 to a -H). This below is for when I was charging things from a -H uncharged to -H charged, although it also happens when I'm actually converting the -CH3 to -H (at Lambdas greater than 85%). I made sure to keep the charges balanced at 0 while mutating and I did it at 3 steps like Michael Shirts suggested. Set 1: turn R-CH3 charges off in a way that preserves the total charge. Set 2: change CH3 LJ to H LJ Set 3: Turn R-H charges on in a way that preserves the total charge. In the portion of the error below atom 9 is -C9-(H67,H68,H66) which in this specific case H67 is already a dummy molecule with no mass or charge. From what I can see, it seems that the atoms do not know how to treat the dummy molecules in terms of angles. How should I treat the dummy molecules? Should I be treating them like hollow spheres with no charge so I would assign them angle constraints? I think it can also be that I'm peturbing the dihedral angles incorrectly. I received errors at first saying that dihedral multiplicities can't be peturbed so I had to equal the multiplicities just to get it to run. Does anybody have any experience with this? Thank you for your help. -Fabian Casteblanco Portion of Error Output: - Reading file nvt0.5.tpr, VERSION 4.5.3 (single precision) starting mdrun 'SIMVASTATIN' 15 steps,300.0 ps. Step 7, time 0.014 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.006111, max 0.139443 (between atoms 9 and 67) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length Step 8, time 0.016 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.007341, max 0.167622 (between atoms 9 and 67) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length Step 8, time 0.016 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.008771, max 0.201182 (between atoms 9 and 67) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length On Mon, Oct 10, 2011 at 4:10 PM, Fabian Casteblanco wrote: > Hi all, > > I have an additional question related to what Steven Neumann was > mentioning. I actually have to do a molecule mutation. I'm trying to > use Michael Shirts method 1) making small > changes 'alchemical' changes in the molecules and computing the free > energies by any method (BAR, TI, etc). I'm specifically want to try > to use BAR at the end once I collect all the data. This helped a lot > on clarification since it seemed that Justin's tutorial is essentially > a FEP except its using the BAR mathematical method for computing the > complete decoupling of the molecule rather than using the old FEP > mathematics of the exponential averaging formula. So BAR is only > referring to the mathematical code used to calculate the overall free > energy for the FEP, correct? > > My question is, for a mutation of a -CH3 group to a -H group, is it > better to simply run: > [+ from (Lambda=0 , R-CH3, full charges and interactions -STATE A) > --> (Lambda=1, R-CH, full charges and interactions -STATE B)] > > OR > > [1) from (Lambda=0 , R-CH3, STATE A : Charges and LJ Interactions: ON) > 2) (Lambda=?, -CH3 Charges: OFF ,LJ Interactions: ON and unmutated) > 3) (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ Interactions: OFF) > 4) (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ interactions: ON and Mutated to > -H) > 5) (Lambda=1, R-CH3, -CH3 STATE B : Charges and LJ Interactions: ON) > > Reason I'm asking is because when I try the first choice to do it > STATE A to STATE B in one step, when I reach Lambda=0.85 and above on > the NVT equilibration right after EM, I receive errors saying that > bonds are moving way to far off their constraints which leads me to > believe that the system is moving too far from where it was energy > minimized. Errors such as: > > Step 188, time 0.376 (ps) > LINCS WARNING > relative constraint deviation after LINCS: > rms 0.17, max 0.000636 (between atoms 9 and 68) > bonds that rotated more than 30 degrees: > atom 1 atom 2 angle previous, current, constraint length > 9 68 31.2 0. 0.1110 0. > > Step 188, time 0.376 (ps) LINCS WARNING > relative constraint deviation after LINCS: > rms 0.15, max 0.000531 (between atoms 9 and 68) > bonds that rotated more than 30 degrees: > atom 1 atom 2 angle previous, current, constraint length > 9 68 31.0 0. 0.1110 0. > > > **Please, if anybody can help, I would greatly appreciate it. Thanks. > -- > Best regards, > > F
[gmx-users] Re: FEP
Thanks Michael for your help. This really helps a lot. Thank you! On Mon, Oct 10, 2011 at 4:10 PM, Fabian Casteblanco wrote: > Hi all, > > I have an additional question related to what Steven Neumann was > mentioning. I actually have to do a molecule mutation. I'm trying to > use Michael Shirts method 1) making small > changes 'alchemical' changes in the molecules and computing the free > energies by any method (BAR, TI, etc). I'm specifically want to try > to use BAR at the end once I collect all the data. This helped a lot > on clarification since it seemed that Justin's tutorial is essentially > a FEP except its using the BAR mathematical method for computing the > complete decoupling of the molecule rather than using the old FEP > mathematics of the exponential averaging formula. So BAR is only > referring to the mathematical code used to calculate the overall free > energy for the FEP, correct? > > My question is, for a mutation of a -CH3 group to a -H group, is it > better to simply run: > [+ from (Lambda=0 , R-CH3, full charges and interactions -STATE A) > --> (Lambda=1, R-CH, full charges and interactions -STATE B)] > > OR > > [1) from (Lambda=0 , R-CH3, STATE A : Charges and LJ Interactions: ON) > 2) (Lambda=?, -CH3 Charges: OFF ,LJ Interactions: ON and unmutated) > 3) (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ Interactions: OFF) > 4) (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ interactions: ON and Mutated to > -H) > 5) (Lambda=1, R-CH3, -CH3 STATE B : Charges and LJ Interactions: ON) > > Reason I'm asking is because when I try the first choice to do it > STATE A to STATE B in one step, when I reach Lambda=0.85 and above on > the NVT equilibration right after EM, I receive errors saying that > bonds are moving way to far off their constraints which leads me to > believe that the system is moving too far from where it was energy > minimized. Errors such as: > > Step 188, time 0.376 (ps) > LINCS WARNING > relative constraint deviation after LINCS: > rms 0.17, max 0.000636 (between atoms 9 and 68) > bonds that rotated more than 30 degrees: > atom 1 atom 2 angle previous, current, constraint length > 9 68 31.2 0. 0.1110 0. > > Step 188, time 0.376 (ps) LINCS WARNING > relative constraint deviation after LINCS: > rms 0.15, max 0.000531 (between atoms 9 and 68) > bonds that rotated more than 30 degrees: > atom 1 atom 2 angle previous, current, constraint length > 9 68 31.0 0. 0.1110 0. > > > **Please, if anybody can help, I would greatly appreciate it. Thanks. > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering > C: +908 917 0723 > E: fabian.castebla...@gmail.com > -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] FEP
Hi all, I have an additional question related to what Steven Neumann was mentioning. I actually have to do a molecule mutation. I'm trying to use Michael Shirts method 1) making small changes 'alchemical' changes in the molecules and computing the free energies by any method (BAR, TI, etc). I'm specifically want to try to use BAR at the end once I collect all the data. This helped a lot on clarification since it seemed that Justin's tutorial is essentially a FEP except its using the BAR mathematical method for computing the complete decoupling of the molecule rather than using the old FEP mathematics of the exponential averaging formula. So BAR is only referring to the mathematical code used to calculate the overall free energy for the FEP, correct? My question is, for a mutation of a -CH3 group to a -H group, is it better to simply run: [+ from (Lambda=0 , R-CH3, full charges and interactions -STATE A) --> (Lambda=1, R-CH, full charges and interactions -STATE B)] OR [1) from (Lambda=0 , R-CH3, STATE A : Charges and LJ Interactions: ON) 2) (Lambda=?, -CH3 Charges: OFF ,LJ Interactions: ON and unmutated) 3) (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ Interactions: OFF) 4) (Lambda=?, R-CH3, -CH3 Charges: OFF ,LJ interactions: ON and Mutated to -H) 5) (Lambda=1, R-CH3, -CH3 STATE B : Charges and LJ Interactions: ON) Reason I'm asking is because when I try the first choice to do it STATE A to STATE B in one step, when I reach Lambda=0.85 and above on the NVT equilibration right after EM, I receive errors saying that bonds are moving way to far off their constraints which leads me to believe that the system is moving too far from where it was energy minimized. Errors such as: Step 188, time 0.376 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.17, max 0.000636 (between atoms 9 and 68) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length 9 68 31.20. 0.1110 0. Step 188, time 0.376 (ps) LINCS WARNING relative constraint deviation after LINCS: rms 0.15, max 0.000531 (between atoms 9 and 68) bonds that rotated more than 30 degrees: atom 1 atom 2 angle previous, current, constraint length 9 68 31.00. 0.1110 0. **Please, if anybody can help, I would greatly appreciate it. Thanks. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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: Question about Justin's Free Energy Tutorial
Thank you guys for your help. I appreciate it. On Thu, Oct 6, 2011 at 1:05 PM, Fabian Casteblanco wrote: > Hello Justin, > > I have a question about your tutorial. If I want to mutate one small > group of a molecule, I would have to not provide 'couple_lambda0' and > 'couple_lambda1', correct? I would essentially have to follow sec > 5.7.4 in the Gromacs manual and I have to actually provide all state A > variable and all state B variables. Gromacs would calculate the new B > state parameters for bond lengths, angles, etc, correct? Are there > any other major differences to account for? > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > E: fabian.castebla...@gmail.com > -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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 Justin's Free Energy Tutorial
Hello Justin, I have a question about your tutorial. If I want to mutate one small group of a molecule, I would have to not provide 'couple_lambda0' and 'couple_lambda1', correct? I would essentially have to follow sec 5.7.4 in the Gromacs manual and I have to actually provide all state A variable and all state B variables. Gromacs would calculate the new B state parameters for bond lengths, angles, etc, correct? Are there any other major differences to account for? -- Best regards, Fabian F. Casteblanco Rutgers University -- E: fabian.castebla...@gmail.com -- 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 Question
Hello all, I have a general question about calculating free energies. I recently used g_bar to calculate the free energies of decoupling coulombic and vdW forces of a solute molecule in solvent. I now need to calculate the free energy of a solute molecule mutating to a new molecule (identical but with an extra -CH3 group) in a vacuum and in a solvent. Does anybody know if it would be better to use g_Bar again but this time having the extra -CH3 group completely decoupled by itself or with it be better to use the older slow-growth method that the Gromacs Manual talks about in Sec 3.12. I'm not sure if there is a way to only decouple a certain part of a molecule using the g_bar method, rather than decoupling the entire molecule. Thanks for your help. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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 sampling using G_bar
Hello all, I am finished running a free energy calculation using g_bar and i followed Justin Lemkuls tutorial and I am in the process of analyzing inorder to determine if I had adequate sampling. I have read the 'BAR' paper by Bennett but there are still some concerns whether I have enough sampling (see attached). I see in the tutorial that the sampling goes up to very large sample numbers, some greater than 30,000 where mine only make it to about 300. This 'samples' refers to how many times a certain energy level was experienced? Does this mean my system should be left for longer times? My lines don't line as as closely as the Methane Decoupling tutorial do. Also, for other solvents that I used, some are giving me a warning for violating the 2nd law of thermodynamics and i noticed its because on one of the output lines, the entropy decreases a small bit, and some are giving me a free energy with a std deviation of at most 1.7 kJ/mol. Does this all mean my sampling is not enough? Thanks for your help. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com H2.agr.tar.gz Description: GNU Zip compressed data -- 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: Free energy calculation question
Thanks Justin! On Wed, Aug 31, 2011 at 11:24 AM, Fabian Casteblanco wrote: > Hello Justin, > > I'm calculating the free energy of a drug in an alcohol solvent. I > have a question referring to your free energy tutorial. You mentioned > that decoupling of electrostatic interactions is linear and decoupling > of vdW can vary. Is this true for your case of methanol in water or > for all cases? When I ran my system of drug in ethanol solvent, I got > a non linear dG for both electrostatic and vdW. Also, is there no > need to find dG of cav ( the free energy required to form the solute > cavity within the solvent) ? I have attached some graphs. > > Thanks for your help. Your tutorial was extremely useful. > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering PhD Student > C: +908 917 0723 > E: fabian.castebla...@gmail.com > -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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: Free energy calculation question
d sc-sigma= 0.3 ;*** - the soft-core sigma for particles which have a C6 or C12 parameter equal to zero or a sigma ;smaller than sc_sigma couple-moltype = LOV ;*** - name of moleculetype to decouple couple-lambda0 = vdw-q ;*** - all interactions are on at lambda=0 couple-lambda1 = vdw ;*** - only vdw interactions are on at lambda=1 couple-intramol = no;*** - Do not decouple intramolecular interactions. That is, the λ factor is applied to only solute-solvent ;nonbonded interactions and not solute-solute nonbonded interactions. nstdhdl = 100 ;*** - The frequency with which ∂H/∂λ and ΔH are written to dhdl.xvg. A value of 100 would probably suffice, since ;the resulting values will be highly correlated and the files will get very large. You may wish to increase this ;value to 100 for your own work. ;Velocity generation gen_vel =no ;Velocity generation is off ;END On Wed, Aug 31, 2011 at 11:24 AM, Fabian Casteblanco wrote: > Hello Justin, > > I'm calculating the free energy of a drug in an alcohol solvent. I > have a question referring to your free energy tutorial. You mentioned > that decoupling of electrostatic interactions is linear and decoupling > of vdW can vary. Is this true for your case of methanol in water or > for all cases? When I ran my system of drug in ethanol solvent, I got > a non linear dG for both electrostatic and vdW. Also, is there no > need to find dG of cav ( the free energy required to form the solute > cavity within the solvent) ? I have attached some graphs. > > Thanks for your help. Your tutorial was extremely useful. > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering PhD Student > C: +908 917 0723 > E: fabian.castebla...@gmail.com > -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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: Free energy calculation question
Hello Justin, I'm looking at the dG vs Lambda plot that is an output from G_bar but on the Shirts et al paper that you included, the part where it talks about linearity, it is referring to dH/dLambda for electrostatic decoupling. So if I take the line formed by dGtotal vs Lambda from g_bar output and take the derivative excluding the beginning and end parts (lambda approaches 0 and 1) then it does seem to take more of a somewhat linear shape. If you can see attached, I think perhaps thats how I was misreading it? I hope that can be the reason, since I've done it with 4 different cases (2 different amount of solvent molecules and using MD and SD integrator) and they all seem to be giving similar results. Thanks for your help. -Fabian On Wed, Aug 31, 2011 at 11:24 AM, Fabian Casteblanco wrote: > Hello Justin, > > I'm calculating the free energy of a drug in an alcohol solvent. I > have a question referring to your free energy tutorial. You mentioned > that decoupling of electrostatic interactions is linear and decoupling > of vdW can vary. Is this true for your case of methanol in water or > for all cases? When I ran my system of drug in ethanol solvent, I got > a non linear dG for both electrostatic and vdW. Also, is there no > need to find dG of cav ( the free energy required to form the solute > cavity within the solvent) ? I have attached some graphs. > > Thanks for your help. Your tutorial was extremely useful. > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering PhD Student > C: +908 917 0723 > E: fabian.castebla...@gmail.com > -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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 calculation question
Hello Justin, I'm calculating the free energy of a drug in an alcohol solvent. I have a question referring to your free energy tutorial. You mentioned that decoupling of electrostatic interactions is linear and decoupling of vdW can vary. Is this true for your case of methanol in water or for all cases? When I ran my system of drug in ethanol solvent, I got a non linear dG for both electrostatic and vdW. Also, is there no need to find dG of cav ( the free energy required to form the solute cavity within the solvent) ? I have attached some graphs. Thanks for your help. Your tutorial was extremely useful. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com dG300_summary.agr Description: application/grace dG300_summary.agr Description: application/grace -- 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: Convert drug Charmm topology to Gromacs
Hello Dr. Alexandre Suman de Araujo, I am not sure how to upload to Gromacs User Contribution. I see downloadable files but not place to upload. Do you have the link where I can upload the files? Thanks, Fabian On Wed, Aug 24, 2011 at 2:55 PM, Fabian Casteblanco wrote: > Hello Steven Neumann, > > I recently converted CGenFF parameters into files that are used by > Gromacs. If this is what you need, shoot me an email and I can > provide you with the data sets. > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering PhD Student > C: +908 917 0723 > E: fabian.castebla...@gmail.com > -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Convert drug Charmm topology to Gromacs
Hello Steven Neumann, I recently converted CGenFF parameters into files that are used by Gromacs. If this is what you need, shoot me an email and I can provide you with the data sets. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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
Dear all, I was running free energy calculation for a drug molecule in solvent. First, For [coulomb + vdW] --> [vdW] , I used 'md' integrator For [vdW] --> [none], I was using 'md' but it required me to switch to 'sd' based on this error message: "WARNING: For proper sampling of the (nearly) decoupled state, stochastic dynamics should be used. " Is using sd really much more accurate? Is it best to use 'sd' for [coulomb + vdW] --> [vdW] as well? Thank you. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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 size selection
Hello all, I'm running a free energy simulation using BAR method for a one drug molecule in several alcohols. I started off using 500 solvent molecules for but realized it was taking too long for the smaller solvents, which would make the larger solvents even more time consuming. I reduced it to two test cases. One using 210 solvent molecules and another using 300 solvent molecules. From what I've read on the user archives, the pressure will always oscillate and for smaller systems, the oscillations will be even larger. For example, on my 210 solvent case, after 1000 ps, the pressure is 1.7 bar (reference set to 1 bar). Is this due to the pressure oscillations? Will running for even longer times eventually make the error at least a bit smaller? Is this ok for free energy calculations? Thanks for your expertise. Your help is appreciated. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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: Free energy calculation
Thanks Justin. Best regards, Fabian On Mon, Aug 8, 2011 at 5:49 PM, Fabian Casteblanco wrote: > Hello all, > > I am setting up a free energy calculation (drug from full coulomb+vdW > in solution --> drug with only vdW in solution --> dummy drug in > solution). > > After reading most of the papers, I understand that you need > significant overlap from the energies for each intermediate point to > overlap so its best to have many intermediate points from Lambda=0 to > Lambda=1. The drug molecule I have is a bit complex but I wasn't sure > if using too small of an intermediate could have a bad effect on the > free_energy calculation. I know Justin Lemkul said in its tutorial > that Lambda=+0.05 should be good for most but I decided to go with > Lambda=+0.02. Could this have any negative effect other than taking a > longer time? Also, how does one come up with the best soft-core > potential parameters to use? Is it ok to use the one from the Methane > in water tutorial? > > Thanks for everyone's expertise. > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering PhD Student > C: +908 917 0723 > E: fabian.castebla...@gmail.com > -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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 calculation
Hello all, I am setting up a free energy calculation (drug from full coulomb+vdW in solution --> drug with only vdW in solution --> dummy drug in solution). After reading most of the papers, I understand that you need significant overlap from the energies for each intermediate point to overlap so its best to have many intermediate points from Lambda=0 to Lambda=1. The drug molecule I have is a bit complex but I wasn't sure if using too small of an intermediate could have a bad effect on the free_energy calculation. I know Justin Lemkul said in its tutorial that Lambda=+0.05 should be good for most but I decided to go with Lambda=+0.02. Could this have any negative effect other than taking a longer time? Also, how does one come up with the best soft-core potential parameters to use? Is it ok to use the one from the Methane in water tutorial? Thanks for everyone's expertise. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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: Dihedral Question
Thanks Justin for your help :-) Best regards, Fabian Casteblanco On Fri, Aug 5, 2011 at 11:46 AM, Fabian Casteblanco wrote: > Hello all, > > I'm having trouble finding information that explains the relationship > between 1,4 pair interactions and dihedrals. I'm building a drug > molecule using CGenFF, and extension of CHARMM27 and since the drug is > a bit complex, it has several dihedral parameters that are missing. I > ran them with the missing parameters and then I reran it with > estimated parameters based on other dihedral parameters for similar > dihedrals and they look somewhat similar but still have noticeable > differences. > > From what I understood, the 1,4 pair interactions should already be > setting most of the dihedral angles up and the dihedrals are simply > added energy consequences inorder to keep the specified torsion angle. > Is this the correct way 1,4 pair interactions and dihedrals work > together? Also, if I were to simply run without using the dihedral > parameters, would it make a huge difference when extracting > themodynamics information from a mixture of this single drug molecule > in different types of solvent? > > Thank you all for sharing your expertise. It's greatly appreciated. > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering PhD Student > C: +908 917 0723 > E: fabian.castebla...@gmail.com > -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Dihedral Question
Hello all, I'm having trouble finding information that explains the relationship between 1,4 pair interactions and dihedrals. I'm building a drug molecule using CGenFF, and extension of CHARMM27 and since the drug is a bit complex, it has several dihedral parameters that are missing. I ran them with the missing parameters and then I reran it with estimated parameters based on other dihedral parameters for similar dihedrals and they look somewhat similar but still have noticeable differences. >From what I understood, the 1,4 pair interactions should already be setting most of the dihedral angles up and the dihedrals are simply added energy consequences inorder to keep the specified torsion angle. Is this the correct way 1,4 pair interactions and dihedrals work together? Also, if I were to simply run without using the dihedral parameters, would it make a huge difference when extracting themodynamics information from a mixture of this single drug molecule in different types of solvent? Thank you all for sharing your expertise. It's greatly appreciated. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] vdW cutoff
Hello, I am quite confused on whether it is better to use a standard cut-off scheme for vdW interactions or if its better to use a switch or shift function for this. I am doing a free energy calculation on the solvation of a drug molecule in a solvent (on CHARMM ff) so I want to be as accurate as possible. The CHARMM paper does state that they use some type of switch function between 10 A and 12 A but I'm confused on the difference between 'Switch' and 'Shift'. I ran the solvents alone using these different types and it seems to give similar results. Is this a major decision to getting accurate free energy calculations? Will using standard 14 A cutoffs with DispersionCorrection=EnerPres be enough? Thanks for anyone with knowledge in this area. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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: Using new atom types
Hello Justin,
Thanks for your help.
Your response:
No, you need the [defaults] directive from the force field. If you
want to make
a completely custom entity, then build your own, although likely most of the
highest-level force field files are going to be the same. You can always add a
new [atomtypes] directive in a proper location in your topology.
Comment: So you are saying I need to refer to the [defaults]
directive so I need to include {#include "charmm27.ff/forcefield.itp"}
on my topology file? Below is the [defaults] directive for CHARMM. I
can simply add my...
>> -CgenFFbon.itp (bonded, angle, dihedral parameters,etc)
>> -CgenFFnb.itp (nonbonded LJ values)
>> -CgenFF.atp(atom type file)
to this section?
* Charmm to Gromacs port writted by *
* Par Bjelkmar, bjelk...@cbr.su.se *
* Alternative correspondance:*
* Erik Lindahl, lind...@cbr.su.se*
#define _FF_CHARMM
[ defaults ]
; nbfunccomb-rule gen-pairs fudgeLJ fudgeQQ
1 2 yes 1.0 1.0
#include "ffcharmm27nb.itp"
#include "ffcharmm27bon.itp"
**???{ #include add new CgenFFbon.itp and CgenFFnb.itp here? }???**?
Or would I need to create one almost identical to this except
replacing the old ffcharm27.itp with my new Cgen.itp files?
This is my sample topology. Your last statement said I could add
[atomtypes] to this section. So I can simply include the [defaults]
directive but just add the Cgenffbon.itp and cCgenffnb.itp to this
section here?
;SAMPLE topology:
;
;
; Include forcefield parameters
#include "charmm27.ff/forcefield.itp"
**???{ #include add new CgenFFbon.itp and CgenFFnb.itp here? }???**?
#include "drug.itp"
[ system ]
; Name
DrugName
[ molecules ]
; Compound#mols
DDD 1
Thanks again for your help. Greatly appreciated!!
Fabian Casteblanco
Rutgers University
--
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: Using new atom types
Hello Peter C. Lai, Thanks for your help. I read your back and forth posts on the GMX Forums with Mark Abraham. Is there any way that I can simply use it from the topology without having to include it on the charmm27.ff/forcefield.itp? I have to run this on a cluster and I don't think I have access to change internal files. I do not need to run using pdb2gmx. Also, how did you manage to include the LJ 1-4 parameters? Did you have to write out each and every pair on GROMACS? Thanks for you help. -Fabian On Thu, Jul 7, 2011 at 2:16 PM, Fabian Casteblanco wrote: > Hello all, > > I'm building a drug molecule using CHARMM parameters. The thing is > that there is this new CGenFF (an extension of CHARMM, but very > similar to the old CHARMM atom types > http://mackerell.umaryland.edu/~kenno/cgenff/) that uses better > parameters for drug-like molecules. I would like to use these new > parameters for this drug molecule (pretty big molecule). I took all > the parameter values from CGenFF and I converted them to values used > by Gromacs and I formatted it just like the original ffcharm27.itp > files found in Gromacs. The only thing is that when it comes to > Lennard Jones,1-4 parameters, CHARMM gives you the actual individual > value by atom whereas Gromacs has already all the pair listed values > stored. Is there any way that I can place these 1-4 parameters in the > drug.itp file itself and have Gromacs calculate the combinations just > as it does for the regular LJ values? > > I have: > -CgenFFbon.itp (bonded, angle, dihedral parameters,etc) > -CgenFFnb.itp (nonbonded LJ values) > -CgenFF.atp (atom type file) > > What exactly do I need to include on my drug.itp (or drug.top) file > inorder so it can read these parameters from the files above? > > Thanks for anyones help! Greatly appreciated!! > > ;SAMPLE topology: > ; > ; > ; Include forcefield parameters > #include "charmm27.ff/forcefield.itp" > > #include "drug.itp" > > [ system ] > ; Name > DrugName > > [ molecules ] > ; Compound #mols > DDD 1 > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering PhD Student > C: +908 917 0723 > E: fabian.castebla...@gmail.com > -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Using new atom types
Hello all, I'm building a drug molecule using CHARMM parameters. The thing is that there is this new CGenFF (an extension of CHARMM, but very similar to the old CHARMM atom types http://mackerell.umaryland.edu/~kenno/cgenff/) that uses better parameters for drug-like molecules. I would like to use these new parameters for this drug molecule (pretty big molecule). I took all the parameter values from CGenFF and I converted them to values used by Gromacs and I formatted it just like the original ffcharm27.itp files found in Gromacs. The only thing is that when it comes to Lennard Jones,1-4 parameters, CHARMM gives you the actual individual value by atom whereas Gromacs has already all the pair listed values stored. Is there any way that I can place these 1-4 parameters in the drug.itp file itself and have Gromacs calculate the combinations just as it does for the regular LJ values? I have: -CgenFFbon.itp (bonded, angle, dihedral parameters,etc) -CgenFFnb.itp (nonbonded LJ values) -CgenFF.atp(atom type file) What exactly do I need to include on my drug.itp (or drug.top) file inorder so it can read these parameters from the files above? Thanks for anyones help! Greatly appreciated!! ;SAMPLE topology: ; ; ; Include forcefield parameters #include "charmm27.ff/forcefield.itp" #include "drug.itp" [ system ] ; Name DrugName [ molecules ] ; Compound#mols DDD1 -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Heat Capacity Question
Hello, I am trying to compare the heat capacity (at NPT) of 1000 molecules of methanol. I ran it all to equilibrium and used g_energy and -nmol 1000 to calculate the heat capacity. The value I achieved, ~140 J/mol*K is far from what I see is 79.5 J/mol. I have read on some past posts that there was some bug in calculating heat capacity but I'm not sure if that applies. I tried plotting Enthalpy vs Temp inorder to get dH/dT since that should be equal to heat capacity but I'm getting ~140 J/mol*K. I'm not sure if I need to include the number of constraints and how would I do that. I used constr=all-bonds. Is there something I'm doing wrong? when it says 140 J/mol*K, it is per mole of what? I had calculated heat of vaporization using Eg - Epot+RT=Hvap and I got it close to the experimental value (34.09 compared to 35.2 kJ/mol). Thanks for the help. -- *Best regards,* ** *Fabian F. Casteblanco* *Rutgers University -- * *Chemical Engineering PhD Student* *C: +908 917 0723* *E: **fabian.castebla...@gmail.com* -- 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] Liquid/Gas Systems
Hi Justin, You had helped me earlier on calculating the heat of vaporization of methanol and it worked great. I'm just trying hard to understand conceptually what is the difference between simulating a liquid phase and a gas phase in Gromacs. I mean technically if we throw in 1000 molecules of component A, would it only be a gas if we made the box super huge? If we run NPT on a system, is that technically when we are compounding it together to find a liquid density? I'm just a bit confused on the difference. I've been simulating liquids up to this point so when I ran only NVT on a single gas molecule, I was trying to understand why we only run NVT on it (0 pressure so no pressure coupling). Thanks for your expertise. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Enthalpy of Vaporization
Thanks Justin. It worked and it is only off by 3% from experimental value so that is great. Thanks for your help! -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Enthalpy of Vaporization
Thanks Justin. That clarified a lot. I'm having trouble simulating the 1 molecule of gas methanol. I took the original energy minimized molecule and I'm trying to start it off on NVT using the following *.mdp file. I got rid of PBC (ns_type=simple) and I set rlist=infinite. After correcting some errors, "Can not have dispersion correction with pbc=no", (set dispcorr=no), and setting nstlist=0 since Gromacs says simulating without cut-offs is usually faster with nstlist=0, I now keep getting an error saying "Can not have nstlist<=0 with twin-range interactions". Should nstlist be =0 as it says in the Gromacs manual when simulating with no cut-offs for pbc=no? It says I can not have Ewald with pbc=no which makes sense but I don't know what to replace it with. Thanks Justin. I appreciate your help. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com nvt.mdp Description: Binary data -- 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] Enthalpy of Vaporization
Hello Justin, Thank you for your response. I just wanted to make sure I understand what you meant. I'm assuming that you want the one molecule in the gas phase at its boiling point and you are doing this because you want the molecule alone with no intermolecular interactions? (since heat of vaporization is the energy required to break those interactions amoung liquid molecules) Is that why we ignore kinetic energy? Is the liquid alcohol that I already simulated (liquid methanol, 1 bar, 298 K) also suppose to be at the 338 K boiling temperature? Would I simply run the npt equilibrium again but simply change the temperature from 298 to 338 K? Also, you stated the equation: DHvap = - + RT The first Epot(gas) would be the potential energy for 1 molecule but the Epot(liq) would be the potential energy for 1 mole? Thanks for your help Justin. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Enthalpy of vaporization
Hello all, I am trying to find the enthalpy of vaporization for 7 types of alcohols to try to compare to experimental values. I have all of them simulated to equilibrium and I can use g_energy to view the Total Energy (both kinetic and potential). In order for me to find the enthalpy of vaporization, is it as simple as just running the NPT script again and changing the temperature to that of its boiling point and 1 degree K over that? Would it just be the difference between the enthalpy at boiling pt and 1 degree past boiling? If anyone could possibly lead me in the right direction I would greatly appreciate it. Thanks for your help. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com npt.mdp Description: Binary data md.mdp Description: Binary data -- 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] Genbox question
Hello, I was wondering if someone can help me with a general genbox question. I have been using the command line: genbox -ci octanol.gro -nmol 1000 -box 5 5 5 -o prop_1000.gro to fill a box with 1000 molecules of octanol. With the smaller n-alcohols, it worked fine but as I started using 1-butanol, 1-pentanol, (longer carbon chains) up to octanol, Gromacs sometimes ends the command with a 'Killedin-Solvent overlap' line at the very end. Sometimes increasing the box size helped, but for 1-pentanol, even 7 7 7 box size didn't work. Also, do these killed runs take up some sort of space. My memory space decreases a little bit even on failed runs. Thanks. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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: Unsteady Density
Thanks Justin. I'm currently running the MD run separately so that the average value won't be affected by the unequilibrated values at the beginning of the NPT run where the density started in the 400s and ran up to the ~high 700s. I had done this before when I had nstlist at 5 but I did run NPT for a lot shorter time last time so Im hoping this time it will be better since last time it was still off (0.781 g/cm3 to 0.791 g.cm3). I will also change the compressibility. Do you think increasing the rlist, rcoulomb, rvdw value to about 1.6, 1.8 nm might also work a bit better? Thanks again. On Mon, May 2, 2011 at 2:16 PM, Fabian Casteblanco wrote: > Hello, > > I have been trying to use CHARMM forcefield to simulate 1000 molecules > of Methanol in a box but I seem to fall short on the density every > time I run. At first, I had my rlist at 1 but CHARMM recommends > greater than 1.2-1.4 nm so I changed accordingly, then I changed it so > that the short neighbor list would update more freqently and finally I > even added a greater amount of steps, hoping that the pressure and > density would settle down at some value. > > Below, I had taken and modified from the Lysozyme tutorial. One thing > I noticed is the contraint -algorithm and how the constraints are set > to 'all-bonds'. Is this necessary? After running NVT (leveled off > near 298K, seemed ok), then NPT, for 300,000 steps, the pressure is > close to the reference 1 bar, ~1.1,1.2 bar, but I noticed the density > is still oscillating, not by much, but it still leaves me questioning > why its not settling around to the literature value of 0.791 g/cm3. I > posted the last several lines from the analysis *.xvg file. NPT > script is also below. > > Is there still more steps that I have to run for? The density just > seems to be oscillating too much and although a literature value of > 0.791 g/cm3 is within the data, it doesnt seem to be settling around > that value. I have similar problems for 1-propanol to where my > simulated density is under that of the literature value. Could it be > possible that maybe I should use different NVT, NPT algorithms > (tcoupl, pcoupl)? > > I would really appreciate anyones help. I'm still in my first few > months of learning the program. Thanks! > > 591.80 38039.125000 297.409546 -206.071747 779.428894 > 592.00 38035.421875 294.95 -20.557419 777.816956 > 592.20 37982.812500 298.063324 -463.192047 776.430542 > 592.40 37681.421875 295.430878 36.529854 776.460876 > 592.60 37476.148438 297.193787 -28.961288 775.634644 > 592.80 37611.574219 297.059875 -553.369263 774.096924 > 593.00 37682.277344 288.297974 -293.567963 776.055237 > 593.20 37147.824219 295.985901 -221.989441 781.570679 > 593.40 37404.906250 293.938721 399.093018 788.921570 > 593.60 37339.914062 297.464783 415.891235 793.838440 > 593.80 37909.019531 292.707184 574.818726 796.730286 > 594.00 37369.820312 304.095703 -119.569496 796.235962 > 594.20 37293.445312 294.449005 -63.335049 794.578552 > 594.40 37196.917969 295.323761 104.148239 791.275085 > 594.60 37776.199219 298.004333 28.989655 785.244629 > 594.80 37893.097656 302.397308 -268.065765 779.457886 > 595.00 38123.460938 295.343109 -140.914047 775.682678 > 595.20 38002.054688 302.233124 -91.893478 774.813660 > 595.40 37835.652344 297.808319 441.931885 776.142090 > 595.60 38103.964844 296.369019 490.766754 778.422302 > 595.80 37686.902344 303.358093 -19.726471 781.057922 > 596.00 37586.890625 293.751648 -92.387199 783.724854 > 596.20 37408.414062 300.007385 -288.156036 785.990417 > 596.40 37894.898438 295.720520 346.861206 788.753174 > 596.60 37673.378906 288.311432 89.059952 789.179993 > 596.80 36881.062500 295.265564 70.789131 791.552673 > 597.00 37077.789062 296.261444 72.207787 793.970398 > 597.20 37123.230469 296.511475 -232.434616 793.653564 > 597.40 37426.152344 304.281891 -133.300797 791.791077 > 597.60 37514.332031 305.183197 -259.364136 787.592712 > 597.80 37393.667969 301.326111 -150.162338 785.053711 > 598.00 37263.054688 299.811554 326.851837 785.553894 > 598.20 37203.261719 300.688904 225.526001 788.008667 > 598.40 37511.636719 302.506714 66.299194 790.725281 > 598.60 37460.843750 299.512909 198.943665 793.35 > 598.80 37248.921875 291.575134 -93.283127 793.582703 > 599.00 36870.480469 293.515381 211.069107 792.929443 > 599.20 37398.097656 293.184448 -163.371140 789.150513 > 599.40 37402.
[gmx-users] Re: Unsteady Density
Hello, Thanks for your response. Here are my responses to your points: nstlist=1 is overkill. I will change this back to 5 except for EM like you said. I originially thought that maybe the density not settling was because the neighbor searching frequency had not been updated frequently enough causing molecule density to change. When I performed this specific NPT run, the average was off due to the fact that the density started off in the 400s and rised up to ~high 700s. The last time I had performed this with nstlist=5 and for the actual MD run, I had a density of approximately 0.781 g/cm3 compared to literature value of 0.791 g/cm3 at 298K. I was getting similar results for 1-propanol with the difference much greater. I was hoping that if I could somehow fix the methanol system, I would be able to fix the 1-propanol system. I created the topologies using the OPLSAA *.itp file (to keep the structure correctly, bonds, angles, etc) and I reassigned the values corresponding to CHARMM (attached). At first, I had placed the values straight from the literature found on the CHARMM website, http://mackerell.umaryland.edu/CHARMM_ff_params.html, but I realized they were nearly identical and were infact the same values in the CHARMM.ff files in GROMACS (i.e. ffbonded.itp, ffnonbonded.itp) All I did was reassign the functions to the functions used in the ffbonded.itp and ffnonbonded.itp found in CHARMM files on GROMACS. The charges were all from the CHARMM parameter website that had this molecule already documented (a copy of the lines is placed at the top of the *.itp file attached). I guess the compressibility factor might make a difference. The manual said that the value is fairly close to most liquids so I had simply ignored the difference but I will look into it. After reading more up on it, I was thinking maybe accounting for gradual cutoffs near the ends using rcoulomb_switch and similarly for rvdw. Anything wrong with the way I created the topology? Dr. Vitaly, I also attached the *.itp files like you mentioned. Thanks again for all your help. On Mon, May 2, 2011 at 2:16 PM, Fabian Casteblanco wrote: > Hello, > > I have been trying to use CHARMM forcefield to simulate 1000 molecules > of Methanol in a box but I seem to fall short on the density every > time I run. At first, I had my rlist at 1 but CHARMM recommends > greater than 1.2-1.4 nm so I changed accordingly, then I changed it so > that the short neighbor list would update more freqently and finally I > even added a greater amount of steps, hoping that the pressure and > density would settle down at some value. > > Below, I had taken and modified from the Lysozyme tutorial. One thing > I noticed is the contraint -algorithm and how the constraints are set > to 'all-bonds'. Is this necessary? After running NVT (leveled off > near 298K, seemed ok), then NPT, for 300,000 steps, the pressure is > close to the reference 1 bar, ~1.1,1.2 bar, but I noticed the density > is still oscillating, not by much, but it still leaves me questioning > why its not settling around to the literature value of 0.791 g/cm3. I > posted the last several lines from the analysis *.xvg file. NPT > script is also below. > > Is there still more steps that I have to run for? The density just > seems to be oscillating too much and although a literature value of > 0.791 g/cm3 is within the data, it doesnt seem to be settling around > that value. I have similar problems for 1-propanol to where my > simulated density is under that of the literature value. Could it be > possible that maybe I should use different NVT, NPT algorithms > (tcoupl, pcoupl)? > > I would really appreciate anyones help. I'm still in my first few > months of learning the program. Thanks! > > 591.80 38039.125000 297.409546 -206.071747 779.428894 > 592.00 38035.421875 294.95 -20.557419 777.816956 > 592.20 37982.812500 298.063324 -463.192047 776.430542 > 592.40 37681.421875 295.430878 36.529854 776.460876 > 592.60 37476.148438 297.193787 -28.961288 775.634644 > 592.80 37611.574219 297.059875 -553.369263 774.096924 > 593.00 37682.277344 288.297974 -293.567963 776.055237 > 593.20 37147.824219 295.985901 -221.989441 781.570679 > 593.40 37404.906250 293.938721 399.093018 788.921570 > 593.60 37339.914062 297.464783 415.891235 793.838440 > 593.80 37909.019531 292.707184 574.818726 796.730286 > 594.00 37369.820312 304.095703 -119.569496 796.235962 > 594.20 37293.445312 294.449005 -63.335049 794.578552 > 594.40 37196.917969 295.323761 104.148239 791.275085 > 594.60 37776.199219 298.004333 28.989655 785.244629 > 594.80 37893.097656 302.397308 -268.065765 779.457886 > 595.00 38123.
[gmx-users] Unsteady Density
Hello, I have been trying to use CHARMM forcefield to simulate 1000 molecules of Methanol in a box but I seem to fall short on the density every time I run. At first, I had my rlist at 1 but CHARMM recommends greater than 1.2-1.4 nm so I changed accordingly, then I changed it so that the short neighbor list would update more freqently and finally I even added a greater amount of steps, hoping that the pressure and density would settle down at some value. Below, I had taken and modified from the Lysozyme tutorial. One thing I noticed is the contraint -algorithm and how the constraints are set to 'all-bonds'. Is this necessary? After running NVT (leveled off near 298K, seemed ok), then NPT, for 300,000 steps, the pressure is close to the reference 1 bar, ~1.1,1.2 bar, but I noticed the density is still oscillating, not by much, but it still leaves me questioning why its not settling around to the literature value of 0.791 g/cm3. I posted the last several lines from the analysis *.xvg file. NPT script is also below. Is there still more steps that I have to run for? The density just seems to be oscillating too much and although a literature value of 0.791 g/cm3 is within the data, it doesnt seem to be settling around that value. I have similar problems for 1-propanol to where my simulated density is under that of the literature value. Could it be possible that maybe I should use different NVT, NPT algorithms (tcoupl, pcoupl)? I would really appreciate anyones help. I'm still in my first few months of learning the program. Thanks! 591.80 38039.125000 297.409546 -206.071747 779.428894 592.00 38035.421875 294.95 -20.557419 777.816956 592.20 37982.812500 298.063324 -463.192047 776.430542 592.40 37681.421875 295.430878 36.529854 776.460876 592.60 37476.148438 297.193787 -28.961288 775.634644 592.80 37611.574219 297.059875 -553.369263 774.096924 593.00 37682.277344 288.297974 -293.567963 776.055237 593.20 37147.824219 295.985901 -221.989441 781.570679 593.40 37404.906250 293.938721 399.093018 788.921570 593.60 37339.914062 297.464783 415.891235 793.838440 593.80 37909.019531 292.707184 574.818726 796.730286 594.00 37369.820312 304.095703 -119.569496 796.235962 594.20 37293.445312 294.449005 -63.335049 794.578552 594.40 37196.917969 295.323761 104.148239 791.275085 594.60 37776.199219 298.004333 28.989655 785.244629 594.80 37893.097656 302.397308 -268.065765 779.457886 595.00 38123.460938 295.343109 -140.914047 775.682678 595.20 38002.054688 302.233124 -91.893478 774.813660 595.40 37835.652344 297.808319 441.931885 776.142090 595.60 38103.964844 296.369019 490.766754 778.422302 595.80 37686.902344 303.358093 -19.726471 781.057922 596.00 37586.890625 293.751648 -92.387199 783.724854 596.20 37408.414062 300.007385 -288.156036 785.990417 596.40 37894.898438 295.720520 346.861206 788.753174 596.60 37673.378906 288.311432 89.059952 789.179993 596.80 36881.062500 295.265564 70.789131 791.552673 597.00 37077.789062 296.261444 72.207787 793.970398 597.20 37123.230469 296.511475 -232.434616 793.653564 597.40 37426.152344 304.281891 -133.300797 791.791077 597.60 37514.332031 305.183197 -259.364136 787.592712 597.80 37393.667969 301.326111 -150.162338 785.053711 598.00 37263.054688 299.811554 326.851837 785.553894 598.20 37203.261719 300.688904 225.526001 788.008667 598.40 37511.636719 302.506714 66.299194 790.725281 598.60 37460.843750 299.512909 198.943665 793.35 598.80 37248.921875 291.575134 -93.283127 793.582703 599.00 36870.480469 293.515381 211.069107 792.929443 599.20 37398.097656 293.184448 -163.371140 789.150513 599.40 37402.238281 297.402008 45.650658 786.303711 599.60 37419.816406 298.520508 -147.688004 783.492371 599.80 37497.332031 296.529877 -15.991615 782.294250 600.00 37691.468750 298.839844 -68.580627 781.974182 title =Methanol npt equilibration ;Run parameters integrator =md ;leap-frog algorithm nsteps =30 ;2 * 5 = 100 ps dt =0.002 ;2fs ;Output control nstxout =100;save coordinates every 0.2 ps nstvout =100;save velocities every 0.2 ps nstenergy =100;save energies every 0.2 ps nstlog =100;update log file every 0.2 ps ;Bond parameters continuation =yes ;Restarting after NVT constraint_algorithm=lincs ;holonomic constraints constraints =all-bonds ;all bonds (even heavy atom-H bonds)constraind lincs_iter =1 ;accuracy of LINCS linc
[gmx-users] Re: Potential energy methanol in a box
Thanks again for your response. I appreciate all your help. On Sat, Apr 23, 2011 at 6:11 PM, Fabian Casteblanco wrote: > Hello, > > I'm trying to fill a box with methanol using CHARMM FF Parameters. I > also need to do this for ethanol and 1-propanol. > > For ethanol, each individual molecule had approximately -19 kJ/mol of > potential energy, then I placed 1000 in a box, performed nvt, npt, > etc. End Result: A negative potential energy (it means that its > stable?) and a density of 0.785 g/cm3 compared to actual 0.789 g/cm3 > at 1 atm, 298 K. This result seemed ok. > > For methanol and 1-propanol, each individual molecule had a larger > 24 - 64 kJ/mol respectively, due to Coulomb 1,4 pair interactions (for > example, for methanol between the HA 0.09 and the H 0.43 gives a large > 64 kJ/mol). Again placed 1000 each in 2 separate boxes, performed > nvt,npt,etc on each one. End Result: A still positive potential > energy for each box (shouldn't this be negative after running npt? on > both) and a density that was 0.789 g/cm3 (for 1-propanol)compared to > actual 0.803 g/cm3 1-propanol at 1 atm, 298 K. Methanol density a bit > closer 0.787 compared to 0.791 g/cm3. > > My questions are: > -There should definetly be a positive potential energy after running > npt, correct? Otherwise it wouldn't be stable? I understand that the > 1,4 coulomb interactions should cancel out with nearby molecules. > *For CHARMM, they state that unless otherwise 1,4 parameters exist > already, they use a full scale 1.0 of the regular LJ, coulomb > potential. > > -As for the density, I know it's somewhat close but I'm trying to be > as accurate as possible. I read somewhere that the average cut-off > (rlist) is usually 1.4 for CHARMM and I had mine set to 1. I'm > guessing that having a larger cutoff would contribute more molecules > to the LJ, coulomb which might rise the density a bit to the correct > value? > > I'm still new to Gromacs so I appreciate all of anybodys help. Thanks a > bunch. > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering PhD Student > C: +908 917 0723 > E: fabian.castebla...@gmail.com > -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Potential energy methanol in a box
Hello, I'm trying to fill a box with methanol using CHARMM FF Parameters. I also need to do this for ethanol and 1-propanol. For ethanol, each individual molecule had approximately -19 kJ/mol of potential energy, then I placed 1000 in a box, performed nvt, npt, etc. End Result: A negative potential energy (it means that its stable?) and a density of 0.785 g/cm3 compared to actual 0.789 g/cm3 at 1 atm, 298 K. This result seemed ok. For methanol and 1-propanol, each individual molecule had a larger 24 - 64 kJ/mol respectively, due to Coulomb 1,4 pair interactions (for example, for methanol between the HA 0.09 and the H 0.43 gives a large 64 kJ/mol). Again placed 1000 each in 2 separate boxes, performed nvt,npt,etc on each one. End Result: A still positive potential energy for each box (shouldn't this be negative after running npt? on both) and a density that was 0.789 g/cm3 (for 1-propanol)compared to actual 0.803 g/cm3 1-propanol at 1 atm, 298 K. Methanol density a bit closer 0.787 compared to 0.791 g/cm3. My questions are: -There should definetly be a positive potential energy after running npt, correct? Otherwise it wouldn't be stable? I understand that the 1,4 coulomb interactions should cancel out with nearby molecules. *For CHARMM, they state that unless otherwise 1,4 parameters exist already, they use a full scale 1.0 of the regular LJ, coulomb potential. -As for the density, I know it's somewhat close but I'm trying to be as accurate as possible. I read somewhere that the average cut-off (rlist) is usually 1.4 for CHARMM and I had mine set to 1. I'm guessing that having a larger cutoff would contribute more molecules to the LJ, coulomb which might rise the density a bit to the correct value? I'm still new to Gromacs so I appreciate all of anybodys help. Thanks a bunch. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com minim.mdp Description: Binary data nvt.mdp Description: Binary data npt.mdp Description: Binary data ethanol.itp Description: Binary data methanol.itp Description: Binary data 1propanol.itp Description: Binary data -- 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: Genbox command question
Thank you Justin. Your explanation really helped. I was actually using 4.0.5, not 4.0.1 (my mistake). Your explanations explained everything very well. Thanks again for your help! On Fri, Apr 15, 2011 at 3:33 PM, Fabian Casteblanco wrote: > Hello, > > I'm trying to place 1000 molecules of 1-propanol using CHARMM FF > parameters in a -box 6 6 6. After minimization, I realized that for > a single molecule, the potential was slightly above 0, '+24 kJ/mol' to > be exact, with electrostatic coloumb potential the greatest > contributor. I used the same thing with OPLSA FF parameters and got a > '+3 kJ/mol'. I realize that when I place 1000 of these molecules, the > potential will be very high since the potentials are positive to begin > with. When I placed 1000 molecules using: > > - "genbox -ci propanol.gro -nmol 1000 -box 6 6 6 -o prop1000.gro" > > I resulted with a box full of 1000 molecules in that specific box size > on gromacs 4.0.1. I did the same on the latest gromacs 4.5 but it > ended up taking longer and the last line stating 'killed'. I redid it > with a box 7 7 7 and it worked on gromacs 4.5. However, after I > energy minimized both 1000 molecules (box 6 6 6 on gromacs 4.0.1) and > (box 7 7 7 on gromacs 4.5), the first gave me a large negative > potential energy (-6000 kJ/mol) after taking a long time and the > second box gave me a somewhat large 2*10^3 kJ/mol. I am still new to > Gromacs but I'm confused on why a certain box size worked on one > version of gromacs and the other didnt. I assume that the difference > in potential energies is since the 2nd box has way more space for > these molecules therefore they are not as well equilibrated as the > more compact box of size 6 6 6. Why does my box size that worked on > gromacs4.0.1 work and on 4.5 it states killed as the last line and > only works redoing it with a bigger box size. > > Thank you so much for your time and help ahead of time. > > -- > Best regards, > > Fabian F. Casteblanco > Rutgers University -- > Chemical Engineering PhD Student > C: +908 917 0723 > E: fabian.castebla...@gmail.com > -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Genbox command question
Hello, I'm trying to place 1000 molecules of 1-propanol using CHARMM FF parameters in a -box 6 6 6. After minimization, I realized that for a single molecule, the potential was slightly above 0, '+24 kJ/mol' to be exact, with electrostatic coloumb potential the greatest contributor. I used the same thing with OPLSA FF parameters and got a '+3 kJ/mol'. I realize that when I place 1000 of these molecules, the potential will be very high since the potentials are positive to begin with. When I placed 1000 molecules using: - "genbox -ci propanol.gro -nmol 1000 -box 6 6 6 -o prop1000.gro" I resulted with a box full of 1000 molecules in that specific box size on gromacs 4.0.1. I did the same on the latest gromacs 4.5 but it ended up taking longer and the last line stating 'killed'. I redid it with a box 7 7 7 and it worked on gromacs 4.5. However, after I energy minimized both 1000 molecules (box 6 6 6 on gromacs 4.0.1) and (box 7 7 7 on gromacs 4.5), the first gave me a large negative potential energy (-6000 kJ/mol) after taking a long time and the second box gave me a somewhat large 2*10^3 kJ/mol. I am still new to Gromacs but I'm confused on why a certain box size worked on one version of gromacs and the other didnt. I assume that the difference in potential energies is since the 2nd box has way more space for these molecules therefore they are not as well equilibrated as the more compact box of size 6 6 6. Why does my box size that worked on gromacs4.0.1 work and on 4.5 it states killed as the last line and only works redoing it with a bigger box size. Thank you so much for your time and help ahead of time. -- Best regards, Fabian F. Casteblanco Rutgers University -- Chemical Engineering PhD Student C: +908 917 0723 E: fabian.castebla...@gmail.com -- 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] Pressure coupling problem
Thank you Justin and Peter for your responses. I tried extending the time on the npt equilibration. It helped but not much. My final pressure after the MD run was about 1.15 bar compared to ref_p which was set to 1 bar. Peter, I will try to analyze the potential error using RMSD and Drift. Thanks. On Mon, Apr 11, 2011 at 4:55 PM, Fabian Casteblanco < fabian.castebla...@gmail.com> wrote: > Hi, > > I'm still in my first few months of using Gromacs. I started by creating > an *.itp and *.top file for *Ethanol* using CHARMM force field > parameters. I made the molecule and it looked fine, put 1000 molecules in a > box, energy minimized it to a negative potential energy, viewed it on VMD, > again looks fine. When I started running the NVT script, I set it equal to > a ref_T of 298 K. It equilibrated at the temperature. Then I tried using > an NPT script to equilibrate it to a ref_p of 1 bar. This is where I get > the problem. The output shows the density is close to the actual > experimental value of 0.789 g/cm^3. But for some reason, my pressure never > gets an average of 1 bar. It keeps oscillating, which I understand is > normal, but the average is always 1.3 or 1.4 bar (it seems the longer I let > it run, the larger the average pressure; 1.38 for 50,000 steps,dt=0.002 and > 1.45 for 75,000 steps,dt=0.002). I don't understand why the ref_p of 1 bar > is not working when I run this NPT.mdp script file. My simple goal is to > have 1000 molecules of ethanol using CHARMM ff parameters at 25degC and 1 > bar and somewhere near the experimental density. > > I would really appreciate anybody's help! I'm new to this but I'm eager to > keep getting better. > > Thanks. > > *NVT SCRIPT (this works fine and takes me to 298 K)* > File Edit Options Buffers Tools Help > title =CHARMM ETHANOL NVT equilibration > ;define =-DPOSRES ;position restrain the protein > ;Run parameters > integrator =md ;leap-frog algorithm > nsteps =5 ;2 * 5 = 100 ps > dt =0.002 ;2fs > ;Output control > nstxout =100;save coordinates every 0.2 ps > nstvout =100;save velocities every 0.2 ps > nstenergy =100;save energies every 0.2 ps > nstlog =100;update log file every 0.2 ps > ;Bond parameters > continuation=no ;first dynamics run > constraint_algorithm=lincs ;holonomic constraints > constraints =all-bonds ;all bonds (even heavy atom-H > bonds)constraind > lincs_iter =1 ;accuracy of LINCS > lincs_order =4 ;also related to accuracy > ;Neighborhood searching > ns_type =grid ;search neighboring grid cells > nstlist =5 ;10 fs > rlist =1.0;short-range neighborlist cutoff (in nm) > rcoulomb=1.0;short-range electrostatic cutoff (in nm) > rvdw=1.0;short-range van der Waals cutoff (in nm) > ;Electrostatics > coulombtype =PME;Particle Mesh Ewald for long-range > electrostat\ > ;ics > pme_order =4 ;cubic interpolation > fourierspacing =0.16 ;grid spacing for FFT > ;Temperature coupling is on > tcoupl =V-rescale ;modified Berendsen thermostat > tc_grps =SYSTEM ;two coupling groups - more accurate > tau_t =0.1;0.1 ;time constant, in ps > ref_t =298;25 ;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=25 ;temperature for Maxwell distribution > gen_seed=-1 ;generate a random seed > ;END > > *NPT SCRIPT* > File Edit Options Buffers Tools Help > title =Ethanol npt equilibration > ;define =-DPOSRES ;position restrain the protein > ;Run parameters > integrator =md ;leap-frog algorithm > nsteps =5 ;2 * 5 = 100 ps > dt =0.002 ;2fs > ;Output control > nstxout =100;save coordinates every 0.2 ps > nstvout =100;save velocities every 0.2 ps > nstenergy =100;save energies every 0.2 ps > nstlog =100;update log file eve
[gmx-users] Pressure coupling problem
Hi, I'm still in my first few months of using Gromacs. I started by creating an *.itp and *.top file for Ethanol using CHARMM force field parameters. I made the molecule and it looked fine, put 1000 molecules in a box, energy minimized it to a negative potential energy, viewed it on VMD, again looks fine. When I started running the NVT script, I set it equal to a ref_T of 298 K. It equilibrated at the temperature. Then I tried using an NPT script to equilibrate it to a ref_p of 1 bar. This is where I get the problem. The output shows the density is close to the actual experimental value of 0.789 g/cm^3. But for some reason, my pressure never gets an average of 1 bar. It keeps oscillating, which I understand is normal, but the average is always 1.3 or 1.4 bar (it seems the longer I let it run, the larger the average pressure; 1.38 for 50,000 steps,dt=0.002 and 1.45 for 75,000 steps,dt=0.002). I don't understand why the ref_p of 1 bar is not working when I run this NPT.mdp script file. My simple goal is to have 1000 molecules of ethanol using CHARMM ff parameters at 25degC and 1 bar and somewhere near the experimental density. I would really appreciate anybody's help! I'm new to this but I'm eager to keep getting better. Thanks. NVT SCRIPT (this works fine and takes me to 298 K) File Edit Options Buffers Tools Help title =CHARMM ETHANOL NVT equilibration ;define =-DPOSRES ;position restrain the protein ;Run parameters integrator =md ;leap-frog algorithm nsteps =5 ;2 * 5 = 100 ps dt =0.002 ;2fs ;Output control nstxout =100 ;save coordinates every 0.2 ps nstvout =100 ;save velocities every 0.2 ps nstenergy =100 ;save energies every 0.2 ps nstlog =100 ;update log file every 0.2 ps ;Bond parameters continuation =no ;first dynamics run constraint_algorithm=lincs ;holonomic constraints constraints =all-bonds ;all bonds (even heavy atom-H bonds)constraind lincs_iter =1 ;accuracy of LINCS lincs_order =4 ;also related to accuracy ;Neighborhood searching ns_type =grid ;search neighboring grid cells nstlist =5 ;10 fs rlist =1.0 ;short-range neighborlist cutoff (in nm) rcoulomb =1.0 ;short-range electrostatic cutoff (in nm) rvdw =1.0 ;short-range van der Waals cutoff (in nm) ;Electrostatics coulombtype =PME ;Particle Mesh Ewald for long-range electrostat\ ;ics pme_order =4 ;cubic interpolation fourierspacing =0.16 ;grid spacing for FFT ;Temperature coupling is on tcoupl =V-rescale ;modified Berendsen thermostat tc_grps =SYSTEM ;two coupling groups - more accurate tau_t =0.1 ;0.1 ;time constant, in ps ref_t =298 ;25 ;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 =25 ;temperature for Maxwell distribution gen_seed =-1 ;generate a random seed ;END NPT SCRIPT File Edit Options Buffers Tools Help title =Ethanol npt equilibration ;define =-DPOSRES ;position restrain the protein ;Run parameters integrator =md ;leap-frog algorithm nsteps =5 ;2 * 5 = 100 ps dt =0.002 ;2fs ;Output control nstxout =100 ;save coordinates every 0.2 ps nstvout =100 ;save velocities every 0.2 ps nstenergy =100 ;save energies every 0.2 ps nstlog =100 ;update log file every 0.2 ps ;Bond parameters continuation =yes ;Restarting after NVT constraint_algorithm=lincs ;holonomic constraints constraints =all-bonds ;all bonds (even heavy atom-H bonds)constraind lincs_iter =1 ;accuracy of LINCS lincs_order =4 ;also related to accuracy ;Neighborhood searching ns_type =grid ;search neighboring grid cells nstlist =5 ;10 fs rlist =1.0 ;short-range neighborlist cutoff (in nm) rcoulomb =1.0 ;short-range electrostatic cutoff (in nm) rvdw =1.0 ;short-range van der Waals cutoff (in nm) ;Electrostatics coulombtype =PME ;Particle Mesh Ewald for long-range electrostat\ ;ics pme_order =4 ;cubic interpolation fourierspacing =0.16 ;grid spacing for FFT ;Temperatur
[gmx-users] Pressure coupling problem
Hi, I'm still in my first few months of using Gromacs. I started by creating an *.itp and *.top file for *Ethanol* using CHARMM force field parameters. I made the molecule and it looked fine, put 1000 molecules in a box, energy minimized it to a negative potential energy, viewed it on VMD, again looks fine. When I started running the NVT script, I set it equal to a ref_T of 298 K. It equilibrated at the temperature. Then I tried using an NPT script to equilibrate it to a ref_p of 1 bar. This is where I get the problem. The output shows the density is close to the actual experimental value of 0.789 g/cm^3. But for some reason, my pressure never gets an average of 1 bar. It keeps oscillating, which I understand is normal, but the average is always 1.3 or 1.4 bar (it seems the longer I let it run, the larger the average pressure; 1.38 for 50,000 steps,dt=0.002 and 1.45 for 75,000 steps,dt=0.002). I don't understand why the ref_p of 1 bar is not working when I run this NPT.mdp script file. My simple goal is to have 1000 molecules of ethanol using CHARMM ff parameters at 25degC and 1 bar and somewhere near the experimental density. I would really appreciate anybody's help! I'm new to this but I'm eager to keep getting better. Thanks. *NVT SCRIPT (this works fine and takes me to 298 K)* File Edit Options Buffers Tools Help title =CHARMM ETHANOL NVT equilibration ;define =-DPOSRES ;position restrain the protein ;Run parameters integrator =md ;leap-frog algorithm nsteps =5 ;2 * 5 = 100 ps dt =0.002 ;2fs ;Output control nstxout =100;save coordinates every 0.2 ps nstvout =100;save velocities every 0.2 ps nstenergy =100;save energies every 0.2 ps nstlog =100;update log file every 0.2 ps ;Bond parameters continuation=no ;first dynamics run constraint_algorithm=lincs ;holonomic constraints constraints =all-bonds ;all bonds (even heavy atom-H bonds)constraind lincs_iter =1 ;accuracy of LINCS lincs_order =4 ;also related to accuracy ;Neighborhood searching ns_type =grid ;search neighboring grid cells nstlist =5 ;10 fs rlist =1.0;short-range neighborlist cutoff (in nm) rcoulomb=1.0;short-range electrostatic cutoff (in nm) rvdw=1.0;short-range van der Waals cutoff (in nm) ;Electrostatics coulombtype =PME;Particle Mesh Ewald for long-range electrostat\ ;ics pme_order =4 ;cubic interpolation fourierspacing =0.16 ;grid spacing for FFT ;Temperature coupling is on tcoupl =V-rescale ;modified Berendsen thermostat tc_grps =SYSTEM ;two coupling groups - more accurate tau_t =0.1;0.1 ;time constant, in ps ref_t =298;25 ;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=25 ;temperature for Maxwell distribution gen_seed=-1 ;generate a random seed ;END *NPT SCRIPT* File Edit Options Buffers Tools Help title =Ethanol npt equilibration ;define =-DPOSRES ;position restrain the protein ;Run parameters integrator =md ;leap-frog algorithm nsteps =5 ;2 * 5 = 100 ps dt =0.002 ;2fs ;Output control nstxout =100;save coordinates every 0.2 ps nstvout =100;save velocities every 0.2 ps nstenergy =100;save energies every 0.2 ps nstlog =100;update log file every 0.2 ps ;Bond parameters continuation=yes;Restarting after NVT constraint_algorithm=lincs ;holonomic constraints constraints =all-bonds ;all bonds (even heavy atom-H bonds)constraind lincs_iter =1 ;accuracy of LINCS lincs_order =4 ;also related to accuracy ;Neighborhood searching ns_type =grid ;search neighboring grid cells nstlist =5 ;10 fs rlist =1.0;short-range neighborlist cutoff (in nm) rcoulomb=1.0;short-range electrostatic cutoff (in nm) rvdw=1.0;short-range van der Waals cutoff (in nm) ;Electrostatics coulombtype =PME;Particle Mesh Ewald for long-range electrostat\ ;ics pme_order =4 ;cubic interpolation fourierspacing =0.16 ;grid spacing for FFT ;Temp
