[gmx-users] Conversion of improper angles from amber object file to gromacs
Hi, I am trying to convert forcefield parameters of zwitterionic glycine from an amber object file format to gromacs .rtp (to use it with amber99sb, which lacks zwitterionic amino acids). The amber library file I plan to use can be found here: https://personalpages.manchester.ac.uk/staff/Richard.Bryce/amber/pro/zaa.off This is what i have so far: [ZG] [ atoms ] N N3 -0.408600 1 H1 H 0.300300 2 H2 H 0.300300 3 H3 H 0.300300 4 CA CT -0.026400 5 HA1 H1 0.073800 6 HA2 H1 0.073800 7 C C 0.835500 8 O O2 -0.724500 9 OXT O2 -0.724500 10 [ bonds ] N H1 N H2 N H3 N CA CA HA1 CA HA2 CA C C O C OXT However, I am unsure about the correct conversion of impropers. In the .off file format, those are based on the atom type, while gromacs expects them to be based on atom names. This leaves room for some ambiguity when the same atom type is used multiple types (i.e the COO or NH3 group). The improper parameter section of the .off file looks like this: !entry.zaa_improper.parm.torsions table str type1 str type2 str type3 str type4 int type dbl kp int n dbl p0 str desc "CT" "O2" "C" "O2" 1 10.50 2 3.141594 "" "CT" "N" "C" "O" 1 10.50 2 3.141594 "" Now, for the CT-O2-C-O2 improper, this would map to two possible combinations: [ impropers ] ; atom1 atom2 atom3 atom4 phase(deg) kd(kJ/mol-1rad-2) dn CA O C OXT 180 43.932 2 CA OXT C O 180 43.932 2 Are both entries required in the gromacs [ improper ] section or must I put only one of them? Thanks, Daniel -- Daniel Bauer, M.Sc. TU Darmstadt Computational Biology & Simulation Schnittspahnstr. 2 64287 Darmstadt ba...@cbs.tu-darmstadt.de -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] No compatible domain decomposition
Hello, I am trying to run a protein simulation consisting of ~100.000 atoms (box dimensions: 8.8x8.8x9.4nm) on multiple mpi threads. For the protein, I need to define a few distance restraints. Right now, these are defined with harmonic bonds in the [ bond ] section of the protein: [ bonds ] ... 53 1952 6 2.437 1000.0 1952 3851 6 2.437 1000.0 3851 5750 6 2.437 1000.0 5750 53 6 2.438 1000.0 53 3851 6 3.447 1000.0 1952 5750 6 3.447 1000.0 Unfortunately, when I try to run this system with more than 1 mpi thread, it fails. Apparently, the high distance for the bond does not allow proper domain decomposition (which makes sense because the bond cannot go over multiple domains). What are my options here? Is there no way to run this across multiple computation nodes except using something like plumed? Initializing Domain Decomposition on 32 ranks Dynamic load balancing: locked Minimum cell size due to atom displacement: 1.301 nm Initial maximum distances in bonded interactions: two-body bonded interactions: 3.444 nm, Harmonic Pot., atoms 53 3851 multi-body bonded interactions: 0.424 nm, Proper Dih., atoms 5136 5145 Minimum cell size due to bonded interactions: 3.789 nm Maximum distance for 5 constraints, at 120 deg. angles, all-trans: 0.855 nm Estimated maximum distance required for P-LINCS: 0.855 nm Using 0 separate PME ranks Scaling the initial minimum size with 1/0.8 (option -dds) = 1.25 Optimizing the DD grid for 32 cells with a minimum initial size of 4.736 nm The maximum allowed number of cells is: X 1 Y 1 Z 1 --- Program: gmx mdrun, version 2019.1 Source file: src/gromacs/domdec/domdec.cpp (line 2403) MPI rank: 0 (out of 32) Fatal error: There is no domain decomposition for 32 ranks that is compatible with the given box and a minimum cell size of 4.73566 nm Change the number of ranks or mdrun option -rdd or -dds Look in the log file for details on the domain decomposition For more information and tips for troubleshooting, please check the GROMACS website at http://www.gromacs.org/Documentation/Errors --- Best regards, Daniel -- Daniel Bauer, M.Sc. TU Darmstadt Computational Biology & Simulation Schnittspahnstr. 2 64287 Darmstadt ba...@cbs.tu-darmstadt.de Don't trust atoms, they make up everything. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Berendsen Barostat for production runs
Hi, I know that the usage of Berendsen baraostat for production runs is generally discouraged because it does not give a correct ensemble average. Nevertheless, I found some recent publications where this is done (presumably for performance reasons). Example: 10.1126/science.1254840 What would be a justification that allows to do this? Under what conditions is it "ok" and when is it not? An explanation attempt by me is that Berendsen is fine if we do not try to investigate properties of the system that depend on the ensemble average (i.e ion conduction events). Best wishes, Daniel -- Daniel Bauer, M.Sc. TU Darmstadt Computational Biology & Simulation Schnittspahnstr. 2 64287 Darmstadt ba...@bio.tu-darmstadt.de Don't trust atoms, they make up everything. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] KCl forms crystals at c=0.15M with amber99sb-ildn
Hi, thanks for the response! Now i know what I have to do :) So did I understand you right that the ion parameters that are currently used in amber99sb by GROMACS are not identical with the ones implemented in the current version of amber99sb used by AMBER? If yes, whats the reason for this (might be more of a question for Justin or Mark). Are there any plans to change this in future releases? Best regards, Daniel On 05/08/2018 04:18 PM, Felipe Merino wrote: > Hi, > > I think the problem comes from the ion parameters within the "vanilla" > amber forcefieCld. If I'm not mistaken, the parameters in the gromacs > port are the Aqvist parameters, which are known to make this little > trick. In the past, we have used the Smith/Dang parameters for > simulations nucleic acids (Dang LX (1995) J Am Chem Soc 117: > 6954–6960) although I think amber is using now Tom Cheatham's > parameters (J. Phys. Chem. B, 2008, 112 (30), pp 9020–9041). > > I hope that helps > Best > Felipe > > On 08/05/18 15:54, Daniel Bauer wrote: >> Hello, >> >> >> We are trying to use amber99sb-ildn with gromacs2018.1 for simulations >> of a protein in a salt solution. However, we observed that KCl seems to >> form stable salt crystals even at concentrations as low as 0.15 mol/l at >> T=298K in our system. Our MDP options have been taken from >> 10.5281/zenodo.1010238 and are listed below. It is reproducible with a >> system containing only 6000 TIP3P water molecules and 16 KCl (0.15M), >> generated by: >> >> gmx genconf -f spc216.gro -o waterbox.gro -nbox 3 3 3 >> >> gmx grompp -f genions -c waterbox -o genions >> >> gmx genion -nname CL -pname K -conc 0.15 -s genions -p topol.top >> >> >> Does someone have an explanation what might cause this to happen? Our >> MDP options are listed below and an RDF showing the accumulation is >> attached. >> >> >> Best Regards, >> >> Daniel >> >> >> Command line: >> ; gmx grompp -f md.mdp -c em/em -o md/md.tpr >> >> ; RUN CONTROL PARAMETERS >> integrator = md >> tinit = 0 >> dt = 0.002 >> nsteps = 125 >> init-step = 0 >> simulation-part = 1 >> comm-mode = Linear >> nstcomm = 1 >> comm-grps = System >> >> ; NEIGHBORSEARCHING PARAMETERS >> cutoff-scheme = Verlet >> nstlist = 10 >> ns-type = Grid >> pbc = xyz >> periodic-molecules = no >> verlet-buffer-tolerance = 0.005 >> rlist = 1.0 >> >> ; OPTIONS FOR ELECTROSTATICS AND VDW >> ; Method for doing electrostatics >> coulombtype = pme >> coulomb-modifier = Potential-shift-Verlet >> rcoulomb-switch = 0 >> rcoulomb = 1.0 >> epsilon-r = 1 >> epsilon-rf = 0 >> vdw-type = Cut-off >> vdw-modifier = Potential-shift-Verlet >> ; cut-off lengths >> rvdw-switch = 0 >> rvdw = 1.0 >> DispCorr = EnerPres >> table-extension = 1 >> energygrp-table = >> fourierspacing = 0.12 >> fourier-nx = 0 >> fourier-ny = 0 >> fourier-nz = 0 >> pme-order = 4 >> ewald-rtol = 1e-05 >> ewald-rtol-lj = 0.001 >> lj-pme-comb-rule = Geometric >> ewald-geometry = 3d >> epsilon-surface = 0 >> >> ; OPTIONS FOR WEAK COUPLING ALGORITHMS >> ; Temperature coupling >> tcoupl = V-rescale >> nsttcouple = -1 >> nh-chain-length = 10 >> tc-grps = System >> tau_t = 0.1 >> ref-t = 298.0 >> >> ; pressure coupling >> Pcoupl = Berendsen >> Pcoupltype = isotropic >> nstpcouple = -1 >> tau_p = 4.0 >> compressibility = 4.5e-5 >> ref_p = 1.0 >> refcoord_scaling = all >> >> ; GENERATE VELOCITIES FOR STARTUP RUN >> gen_vel = yes >> gen_temp = 298.0 >> gen_seed = -1 >> >> ; OPTIONS FOR BONDS >> constraints = all-bonds >> constraint_algorithm = lincs >> continuati
[gmx-users] KCl forms crystals at c=0.15M with amber99sb-ildn
Hello, We are trying to use amber99sb-ildn with gromacs2018.1 for simulations of a protein in a salt solution. However, we observed that KCl seems to form stable salt crystals even at concentrations as low as 0.15 mol/l at T=298K in our system. Our MDP options have been taken from 10.5281/zenodo.1010238 and are listed below. It is reproducible with a system containing only 6000 TIP3P water molecules and 16 KCl (0.15M), generated by: gmx genconf -f spc216.gro -o waterbox.gro -nbox 3 3 3 gmx grompp -f genions -c waterbox -o genions gmx genion -nname CL -pname K -conc 0.15 -s genions -p topol.top Does someone have an explanation what might cause this to happen? Our MDP options are listed below and an RDF showing the accumulation is attached. Best Regards, Daniel Command line: ; gmx grompp -f md.mdp -c em/em -o md/md.tpr ; RUN CONTROL PARAMETERS integrator = md tinit = 0 dt = 0.002 nsteps = 125 init-step = 0 simulation-part = 1 comm-mode = Linear nstcomm = 1 comm-grps = System ; NEIGHBORSEARCHING PARAMETERS cutoff-scheme = Verlet nstlist = 10 ns-type = Grid pbc = xyz periodic-molecules = no verlet-buffer-tolerance = 0.005 rlist = 1.0 ; OPTIONS FOR ELECTROSTATICS AND VDW ; Method for doing electrostatics coulombtype = pme coulomb-modifier = Potential-shift-Verlet rcoulomb-switch = 0 rcoulomb = 1.0 epsilon-r = 1 epsilon-rf = 0 vdw-type = Cut-off vdw-modifier = Potential-shift-Verlet ; cut-off lengths rvdw-switch = 0 rvdw = 1.0 DispCorr = EnerPres table-extension = 1 energygrp-table = fourierspacing = 0.12 fourier-nx = 0 fourier-ny = 0 fourier-nz = 0 pme-order = 4 ewald-rtol = 1e-05 ewald-rtol-lj = 0.001 lj-pme-comb-rule = Geometric ewald-geometry = 3d epsilon-surface = 0 ; OPTIONS FOR WEAK COUPLING ALGORITHMS ; Temperature coupling tcoupl = V-rescale nsttcouple = -1 nh-chain-length = 10 tc-grps = System tau_t = 0.1 ref-t = 298.0 ; pressure coupling Pcoupl = Berendsen Pcoupltype = isotropic nstpcouple = -1 tau_p = 4.0 compressibility = 4.5e-5 ref_p = 1.0 refcoord_scaling = all ; GENERATE VELOCITIES FOR STARTUP RUN gen_vel = yes gen_temp = 298.0 gen_seed = -1 ; OPTIONS FOR BONDS constraints = all-bonds constraint_algorithm = lincs continuation = no lincs_order = 4 lincs_iter = 1 lincs-warnangle = 30 morse = no Best Regards, Daniel -- Daniel Bauer, M.Sc. TU Darmstadt Computational Biology & Simulation Schnittspahnstr. 2 64287 Darmstadt ba...@cbs.tu-darmstadt.de Don't trust atoms, they make up everything. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Minimal PCI Bandwidth for Gromacs and Infiniband?
Hello, In our group, we have multiple identical Ryzen 1700x / Nvidia GeForce 1080 GTX computing nodes and think about interconnecting them via InfiniBands. Does anyone have Information on what Bandwidth is required by GROMACS for communication via InfiniBand (MPI + trajectory writing) and how it scales with the number of nodes? The mainboards we are currently using can only run one PCIe slot with 16 lanes. When using both PICe slots (GPU+InfiniBand), they will run in dual x8 mode (thus bandwidth for both GPU and InfiniBand will be reduced to 8 GB/s instead of 16 GB/s). Now we wonder if the reduced bandwidth will hurt GROMACS performance due to bottlenecks in GPU/CPU communication and/or communication via InfiniBand. If this is the case, we might have to upgrade to new mainboards with dual x16 support. Best regards, Daniel -- Daniel Bauer, M.Sc. TU Darmstadt Computational Biology & Simulation Schnittspahnstr. 2 64287 Darmstadt ba...@cbs.tu-darmstadt.de Don't trust atoms, they make up everything. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Gromacs 2018 beta: g_membed broken?
Hello,
Im using g_membed to insert potassium channels into POPC membranes. My
"standard protocol" which worked very well for me on all proteins I used
so far seems to not work with gromacs 2018 (using the exact same input
files).
This is what I am doing:
1) orienting the protein onto the membrane (lambada)
2) Merging the topology files for the membrane and protein
3) grompp: gmx grompp -f g_membed.mdp -n -c input.gro -r input.gro -o
g_membed.tpr --maxwarn 2
4) gmx mdrun -deffnm g_membed -membed embed.dat -mn index.ndx -mp
embedded.top -c embedded.top # selecting protein and membrane group
First thing I noticed is that the number of warnings i have to ignore
during grompp increased by one because of the new PME warning for
charged systems. My protein has a net charge, so I have to choose if i
insert counter ions before or after protein insertion (and I chose to do
it after insertion).
So the warnings I get are:
WARNING 1: You are using Ewald electrostatics in a system with net
charge // I think I can ignore this since I neutralize right after
insertion.
WARNING 2: Can not exclude the lattice Coulomb energy between energy
groups // g_membed does not support verlet
Finally when I run the insertion simulation on the 2018 beta I get the
following error:
step 0: One or more water molecules can not be settled.
Check for bad contacts and/or reduce the timestep if appropriate.
Back Off! I just backed up step0b.pdb to ./#step0b.pdb.1#
Back Off! I just backed up step0c.pdb to ./#step0c.pdb.1#
Wrote pdb files with previous and current coordinates
step 0Warning: Only triclinic boxes with the first vector parallel to
the x-axis and the second vector in the xy-plane are supported.
Box (3x3):
Box[ 0]={ -nan, -nan, -nan}
Box[ 1]={ -nan, -nan, -nan}
Box[ 2]={ -nan, -nan, -nan}
Can not fix pbc.
This seems to happen independent of my system (I tried several of my
membrane patches and proteins that worked previously). All of them
insert into the membrane with 2016 successfully but fail on the 2018
version. So far I tried to increase the number of steps, reduced the
timestep and played arround with xyinit and zinit as well as setting
-DFLEXIBLE in the mdp, but no success.
Any idea whats wrong? Or should I move this issue directly to redmine
since its related to a version update?
This is my embed.dat file:
nxy = 1000
xyinit = 0.1
xyend = 1.0
rad = 0.22
ndiff = 0
asymmetry = no
pieces = 1
maxwarn = 1
And my mdp:
; Run parameters
integrator = md
nsteps = 1000
dt = 0.002
; I have posres for the protein enabled, although this is irrelevant
because the protein group coords are frozen. Same issue if i disable
this and remove posres!
refcoord_scaling = all
;define = -DFLEXIBLE ; does not help
; membed options
energygrps = insertion_group
freezegrps = insertion_group
freezedim = Y Y Y
energygrp_excl = insertion_group insertion_group
; OUTPUT CONTROL OPTIONS
nstxout = 0
nstvout = 0
nstfout = 0
nstlog = 1000
nstenergy = 100
nstxtcout = 1000
; Bond parameters
continuation = no
constraint_algorithm = lincs
lincs_iter = 2
lincs_order = 8
; CHARMM36 params
constraints = h-bonds
cutoff-scheme = Group ; Verlet is not supported by g_membed
vdwtype = Cut-off
vdw-modifier = force-switch
rlist = 1.2
rvdw = 1.2
rvdw-switch = 0.8
coulombtype = PME
rcoulomb = 1.2
DispCorr = no
; Temperature coupling
tcoupl = V-rescale
tc-grps = Other Water_and_ions Protein
tau_t = 0.1 0.1 0.1
ref-t = 298.0 298.0 298.0
; Pressure coupling
pcoupl = Berendsen
pcoupltype = semiisotropic
tau_p = 5.0
compressibility = 4.5e-5 4.5e-5
ref_p = 1.0 1.0
; Periodic boundary conditions
pbc = xyz
; Velocity generation
gen_vel = yes
gen_temp = 298.0
gen_seed = -1
; COM motion removal
; These options remove motion of the protein/bilayer relative to the
solvent/ions
comm-mode = Linear
comm-grps = Other Water_and_ions Protein
Best regards,
Daniel
--
Daniel Bauer, M.Sc.
TU Darmstadt
Computational Biology & Simulation
Schnittspahnstr. 2
64287 Darmstadt
ba...@cbs.tu-darmstadt.de
Don't trust atoms, they make up everything.
--
Gromacs Users mailing list
* Please search the archive at
http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting!
* Can't post? Read http://www.gromacs.org/Support/Mailing_Lists
* For (un)subscribe requests visit
https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a
mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Conversion of charmm36 parameters from namd to gromacs
Unfortunatly It seems like my university has no contract to access this paper so I cant check the SI. The ones I have are from http://dx.doi.org/10.7554/eLife.25844 (adjusted to better match the free energy of solvation in NMA). On 12/04/2017 07:53 PM, Justin Lemkul wrote: > > > On 12/4/17 12:21 PM, Daniel Bauer wrote: >> Hello, >> >> I finally found my error in the conversion. As always, the devil is in >> the detail. I was under the assumption that parameters listed in the >> original forcefield files (toppar) are also sigma and epsilon values. >> However, as you know, the original files list Rmin values (and not >> sigma). >> >> With the correct conversion term: >> >> eps = Rmin/(10*2^(1/6)) >> >> I can now reproduce the conversion of the numbers from the charmm >> implementation to gromacs (and thus know how to apply my NBFIX for this >> value). > > I don't know what NBFIX you might have, but FYI there is a value that > will become official in the next release of the force field. We've > validated it against proteins and nucleic acids. See the SI of > http://dx.doi.org/10.1128/AAC.01572-17 > > -Justin > -- Daniel Bauer, M.Sc. TU Darmstadt Computational Biology & Simulation Schnittspahnstr. 2 64287 Darmstadt ba...@cbs.tu-darmstadt.de Don't trust atoms, they make up everything. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Conversion of charmm36 parameters from namd to gromacs
Hello, I finally found my error in the conversion. As always, the devil is in the detail. I was under the assumption that parameters listed in the original forcefield files (toppar) are also sigma and epsilon values. However, as you know, the original files list Rmin values (and not sigma). With the correct conversion term: eps = Rmin/(10*2^(1/6)) I can now reproduce the conversion of the numbers from the charmm implementation to gromacs (and thus know how to apply my NBFIX for this value). Thanks for your patience and best regards, Daniel On 12/03/2017 10:18 PM, Daniel Bauer wrote: > Hello, > > > I compared the LJ parameters for the interaction Potassium - backbone > carbonyl oxygen for CHARMM36 between the Gromacs and NAMD version of the > forcefield. I found different numbers for the sigma value i cannot > explain to myself: > > > The stock charmm36 implementation has the following LJ parameters for > potassium and oxygen: > > ; atomEmin (kcal/mol)Rmin/2 (A) > > POT0.0871.76375 > > O0.1201.7 > > > Converting the above values to units used in gromacs (kJ/mol and nm) and > applying standard combination rules this should give the following > nonbonded energy parameters: > > [ pairtypes] > > ; sig = (Rmin,i/2 + Rmin,j/2)/10 > > ; eps = sqrt(Emin,i * Emin,j) * 5.184 kcal/kJ > > ; ijfuncsig (nm)eps (kj/mol) > > OPOT10.346 0.4275 > > However, the actual entry for the interaction in GROMACS/charmm36 has > the energy minimum at a much smaller distance (epsilon is on the point > though!) > > ; ij func sig (nm)eps (kJ/mol) > > OPOT10.2820.4275 > > > Can somebody tell me the reason for this huge difference? Is there an > error in my calculation? I am trying to convert an NBFIX applied to this > specific interaction to gromacs. However, I am not sure how to proceed > without doing the reparametrization in gromacs again, because the stock > values differ that much. > > Best regards, > > Daniel > > -- Daniel Bauer, M.Sc. TU Darmstadt Computational Biology & Simulation Schnittspahnstr. 2 64287 Darmstadt ba...@cbs.tu-darmstadt.de Don't trust atoms, they make up everything. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Conversion of charmm36 parameters from namd to gromacs
Hello Justin, Thanks for your response. I was under the assumption that entries in [ pairtypes ] are the ones that override standard LJ interactions derived by pairgen. I now know this is a different section. From what you are saying I come to the conclusion that the POT O interaction in charmm36 for gromacs have no custom potential but are derived via combination rules (as they are in namd). Nevertheless, the numbers still dont seem to match: Since we can ignore 1-4 Interactions in charmm36 as well as gromacs for POT-O, the resulting parameters from the namd version are: sig = 0.3463750 nm; eps=0.42750571533 kJ/mol However, combining the entries for O and POT in the atomtypes section of the gromacs version i get something different: [ atomtypes ] O 815.9990000.000 A 0.302905564168 0.50208 POT1939.0983000.000 A 0.314264522824 0.36401 === sig = 0.308585043 nm; eps=0.42750689 kJ/mol Again, the epsilon value is matching the namd version while sigma is still off by about 0.4 A. If my conversion of the namd parameters and combination rules are correct, this would be a problem. The different position of the energy minimum results in a potential energy difference of around 7 kJ/mol for a distance of 0.29 nm (typical distance for POT-O in the selectivity filter of potassium channels, see attached image). Can you give me a source where I can read up on why the values for sigma and epsilon in gromacs are the way they are? It must be more than a direct conversion of the numbers from kcal to kJ. Or did I again made a mistake or took the wrong numbers? Best regards, Daniel On 12/03/2017 11:23 PM, Justin Lemkul wrote: > > > On 12/3/17 4:18 PM, Daniel Bauer wrote: >> Hello, >> >> >> I compared the LJ parameters for the interaction Potassium - backbone >> carbonyl oxygen for CHARMM36 between the Gromacs and NAMD version of the >> forcefield. I found different numbers for the sigma value i cannot >> explain to myself: >> >> >> The stock charmm36 implementation has the following LJ parameters for >> potassium and oxygen: >> >> ; atomEmin (kcal/mol)Rmin/2 (A) >> >> POT0.0871.76375 >> >> O0.1201.7 >> >> >> Converting the above values to units used in gromacs (kJ/mol and nm) and >> applying standard combination rules this should give the following >> nonbonded energy parameters: >> >> [ pairtypes] >> >> ; sig = (Rmin,i/2 + Rmin,j/2)/10 >> >> ; eps = sqrt(Emin,i * Emin,j) * 5.184 kcal/kJ >> >> ; ijfuncsig (nm)eps (kj/mol) >> >> OPOT10.346 0.4275 >> >> However, the actual entry for the interaction in GROMACS/charmm36 has >> the energy minimum at a much smaller distance (epsilon is on the point >> though!) >> >> ; ij func sig (nm)eps (kJ/mol) >> >> OPOT10.2820.4275 >> >> >> Can somebody tell me the reason for this huge difference? Is there an >> error in my calculation? I am trying to convert an NBFIX applied to this >> specific interaction to gromacs. However, I am not sure how to proceed >> without doing the reparametrization in gromacs again, because the stock >> values differ that much. > > This is actually a value that isn't used. [pairtypes] are for 1-4 > interactions, which can never occur between a K+ ion and a carbonyl O. > The inclusion in the list is just because it's really hard to tell our > conversion program to exclude some massive list of impossible > interactions. So they get generated. > > But to your point, you're missing a critical point of the CHARMM > parameter format. The full NONBonded parameter entry in > par_all36m_prot.prm for carbonyl O is: > > O 0.00 -0.12 1.70 0.00 -0.12 1.40 > > The last fields are epsilon and Rmin/2 when that atom type is used in > a 1-4 interaction. So -0.12/1.7 is for the normal nonbonded > interaction, -0.12/1.4 is for any 1-4 involving O. > > -Justin > -- Daniel Bauer, M.Sc. TU Darmstadt Computational Biology & Simulation Schnittspahnstr. 2 64287 Darmstadt ba...@cbs.tu-darmstadt.de Don't trust atoms, they make up everything. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Conversion of charmm36 parameters from namd to gromacs
Hello, I compared the LJ parameters for the interaction Potassium - backbone carbonyl oxygen for CHARMM36 between the Gromacs and NAMD version of the forcefield. I found different numbers for the sigma value i cannot explain to myself: The stock charmm36 implementation has the following LJ parameters for potassium and oxygen: ; atomEmin (kcal/mol)Rmin/2 (A) POT0.0871.76375 O0.1201.7 Converting the above values to units used in gromacs (kJ/mol and nm) and applying standard combination rules this should give the following nonbonded energy parameters: [ pairtypes] ; sig = (Rmin,i/2 + Rmin,j/2)/10 ; eps = sqrt(Emin,i * Emin,j) * 5.184 kcal/kJ ; ijfuncsig (nm)eps (kj/mol) OPOT10.346 0.4275 However, the actual entry for the interaction in GROMACS/charmm36 has the energy minimum at a much smaller distance (epsilon is on the point though!) ; ij func sig (nm)eps (kJ/mol) OPOT10.2820.4275 Can somebody tell me the reason for this huge difference? Is there an error in my calculation? I am trying to convert an NBFIX applied to this specific interaction to gromacs. However, I am not sure how to proceed without doing the reparametrization in gromacs again, because the stock values differ that much. Best regards, Daniel -- Daniel Bauer, M.Sc. TU Darmstadt Computational Biology & Simulation Schnittspahnstr. 2 64287 Darmstadt ba...@cbs.tu-darmstadt.de Don't trust atoms, they make up everything. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Lipid Simulation Analysis
Hello, First of all you can try to reproduce basic lipid properties as the area per lipid and bilayer thickness. Both can be easily calculated from the position of lipid headgroups. Other useful properties is the diffusion rate of lipids (accessable by calculating the mean square displacement) and the deuterium order. Best regards! On 31.05.2017 15:39, Mr. Zaved Hazarika wrote: > Hello Nikhil Maroli > > Basically I am new to lipid bilayer simulation. As of now I don't have any > particular objective. I am trying to learn lipid bilayer simulation. > Any suggestion would be great. > > How can I check the stability of my simulated bilayer? > > > ___ > D I S C L A I M E R > This e-mail may contain privileged information and is intended solely for > the individual named. If you are not the named addressee you should not > disseminate, distribute or copy this e-mail. Please notify the sender > immediately by e-mail if you have received this e-mail in error and > destroy it from your system. Though considerable effort has been made to > deliver error free e-mail messages but it can not be guaranteed to be secure > or error-free as information could be intercepted, corrupted, lost, > destroyed, > delayed, or may contain viruses. The recipient must verify the integrity of > this e-mail message. > -- Don't trust atoms, they make up everything. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Continuation mdp flag
Thank you! On 05.05.2017 16:05, Mark Abraham wrote: > Hi, > > It really doesn't do much. > http://manual.gromacs.org/documentation/2016/user-guide/mdp-options.html#bonds > covers > it. (Background: people expect that output is written with constraints > satisfied, we write the checkpoint at the same time for convenience and > efficiency, the initial simulation part probably only has constraints > approximately satisfied from the user input to grompp, and before GROMACS 4 > a .tpr file was the only way to both start and continue a run.) Probably > the mdp option should be removed in favour of logic based around whether an > input checkpoint file (which can be constructed to be sane wrt the intended > simulation) was supplied to mdrun. > > Mark > > On Fri, May 5, 2017 at 3:54 PM Daniel Bauer <ba...@cbs.tu-darmstadt.de> > wrote: > >> Hello, >> >> This might be a trivial question for most of you but can someone give me >> more detail about what the 'continuation' flag in the .mdp settings >> does? The manual does not include a lot of information on that. >> >> Best regards, >> Daniel >> -- >> Gromacs Users mailing list >> >> * Please search the archive at >> http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before >> posting! >> >> * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists >> >> * For (un)subscribe requests visit >> https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or >> send a mail to gmx-users-requ...@gromacs.org. >> -- Don't trust atoms, they make up everything. -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Continuation mdp flag
Hello, This might be a trivial question for most of you but can someone give me more detail about what the 'continuation' flag in the .mdp settings does? The manual does not include a lot of information on that. Best regards, Daniel -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] CHARMM36 and energygroup exclusions
On 05/02/2017 07:45 PM, Justin Lemkul wrote: > > Your input configurations is probably broken across PBC, so grompp > does a clunky calculation and thinks things are split really far > apart. It doesn't matter though, because CHARMM36 doesn't use charge > groups. Every atom is its own group so by definition you can't have > any problems with this. > > -Justin This is very likely the case since the input structure is the output .gro from a previous run equilibration simulation. I think I was irritated because http://www.gromacs.org/Documentation/Errors suggests that for broken molecules "the sum of the two largest charge groups will correspond to a value twice the box vector along which the molecule is broken". However, the sum is not exactly twice my box length in x,y or z. Thinking about it, I guess the vector must not neccessary be in the direction of one axis though. Just to be sure, do you think these settings are fine for the given scenario? Thanks, Daniel -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] CHARMM36 and energygroup exclusions
Hello, I need to run a simulation of a membrane system with CHARMM36 and energygrp exclusions. Therefore, I used the mdp settings suggested by the gromacs page but a lower rvdw-switch (0.8 instead of 1.0) to better reproduce DPPC lipid params: constraints = h-bonds cutoff-scheme = Group vdwtype = cutoff vdw-modifier = force-switch rlist = 1.2 rvdw = 1.2 rvdw-switch = 0.8 coulombtype = PME rcoulomb = 1.2 DispCorr = no Since I need to exclude energy groups, I cannot use the Verlet cufoff scheme but have to use group scheme instead. However, It seems like this incompatible with the cutoff and rlist settings since I get a warning: WARNING 2 [file md0.mdp]: The sum of the two largest charge group radii (19.291325) is larger than rlist (1.20) I think this is because rlist is the same size as rcoulomb and rvdw. However, setting a larger rlist is also not possible because of the PME coulombtype (Error: With coulombtype = PME, rcoulomb must be equal to rlist). My question is is it safe to ignore this warning if the simulation is only runs for a very limited time (2ns) or not? And if not, what changes to my settings would you suggest? To give you some background: I try to add a VdW sphere in my system that only interacts with the protein but allows water to pass through it freely (hence settin energygrp-excl between that particle and everything but the protein). After 2 ns, the sphere will be removed and the system can be simulated with the Verlet scheme and original settings again. Best regards, Daniel B -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] CHARMM36 and energygroup exclusions
Hello, I need to run a simulation of a membrane system with CHARMM36 and energygrp exclusions. Therefore, I used the mdp settings suggested by the gromacs page but a lower rvdw-switch (0.8 instead of 1.0) to better reproduce DPPC lipid params: constraints = h-bonds cutoff-scheme = Group vdwtype = cutoff vdw-modifier = force-switch rlist = 1.2 rvdw = 1.2 rvdw-switch = 0.8 coulombtype = PME rcoulomb = 1.2 DispCorr = no Since I need to exclude energy groups, I cannot use the Verlet cufoff scheme but have to use group scheme instead. However, It seems like this incompatible with the cutoff and rlist settings since I get a warning: WARNING 2 [file md0.mdp]: The sum of the two largest charge group radii (19.291325) is larger than rlist (1.20) I think this is because rlist is the same size as rcoulomb and rvdw. However, setting a larger rlist is also not possible because of the PME coulombtype (Error: With coulombtype = PME, rcoulomb must be equal to rlist). My question is is it safe to ignore this warning if the simulation is only runs for a very limited time (2ns) or not? And if not, what changes to my settings would you suggest? To give you some background: I try to add a VdW sphere in my system that only interacts with the protein but allows water to pass through it freely (hence settin energygrp-excl between that particle and everything but the protein). After 2 ns, the sphere will be removed and the system can be simulated with the Verlet scheme and original settings again. Best regards, Daniel B -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
Re: [gmx-users] Exploding temp/pressure.
Did you forget to set an initial temperature distribution? In your mdp you set gen_vel to no. So it will not create an initial velocity distribution matching your temperature. For the first simulation of your workflow, you should set gen_vel: yes and continuation: no. Then for all simulation that build on top of that simulation you use the velocities from the last timestep of the previous one (mdrun -t flag) and disable gen_vel + set continuation to yes. Also consider running your first simulation with berenden barostat to quickly equilibrate the pressure. For actual sampling of the configuration space you should then switch to parrinello since berendson doesnt sample the ensemble correctly. Both temperature and pressure should equilibrate very quickly in a martini simulation, allowing to increase the timestep to 3 fs or more rather quickly. Hope this helps! Daniel -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
[gmx-users] Should I restrain omega dihedrals or just phi/psi?
Hello, Im working on MD simulations of membrane proteins and want to restrain the backbone dihedrals of some amino acids. I now wonder if its sufficient to restrain only phi/psi dihedrals as shown in the manual ( http://www.gromacs.org/Documentation/How-tos/Dihedral_Restraints), or if I should also set restrains on omega dihedrals? (I dont think it matters but Im using Gromacs 2016.1) Best regards, Daniel Bauer -- Gromacs Users mailing list * Please search the archive at http://www.gromacs.org/Support/Mailing_Lists/GMX-Users_List before posting! * Can't post? Read http://www.gromacs.org/Support/Mailing_Lists * For (un)subscribe requests visit https://maillist.sys.kth.se/mailman/listinfo/gromacs.org_gmx-users or send a mail to gmx-users-requ...@gromacs.org.
