[gmx-users] Re: Good method to turn off electrostatics
Justin A. Lemkul
jalemkul at vt.edu
Tue Jul 12 19:16:09 CEST 2011
Andrew DeYoung wrote:
> Hi,
>
> Many thanks, Justin, for your time. As a first try, I am going ahead with
> your proposed option (1): make my own local copy of spce.itp and adjust the
> call (the #include) in the topology file. I have zeroed the charges in this
> local copy of spce.itp and tried to prepare for energy minimization, using
> grompp:
>
> grompp -f minim.mdp -c water_processed.gro -p water.top -o em.tpr -po
> emout.mdp
>
> When I do this, I get a warning message:
> ------
> WARNING 1 [file minim.mdp]:
> You are using full electrostatics treatment PME for a system without
> charges.
> This costs a lot of performance for just processing zeros, consider using
> Cut-off instead.
> ------
> which leads to a fatal error. I can work around this fatal error by setting
> -maxwarn to 1, and this is what I am going to try now. At this point, I am
> not too concerned about speed because my system is small and my simulation
> is short. But, just out of curiosity, what do you think that the above
> warning is suggesting that I do instead? As it stands now, in my .mdp file,
> I have (among other things) the following:
> ------
> ...
> rlist = 1.0
> coulombtype = PME
> rcoulomb = 1.0
> rvdw = 1.0
> pbc = xyz
> optimize_fft = yes
> ...
> ------
> I guess that it is suggesting that I NOT use coulombtype = PME, and instead
> use coulombtype = Cut-off, as referred to on this page:
Correct.
> http://manual.gromacs.org/current/online/mdp_opt.html#el . But that page
> does seem to say that, even if I use coulombtype = Cut-off, I will still
> need to set rlist and rcoulomb (where rcoulomb >= rlist). My question is, do
> you think, then, that the warning message is suggesting that I use
> coulombtype = Cut-off with rcoulomb = 0.0 and rlist = 0.0?
>
rlist affects neighbor searching and thus van der Waals calculations as well.
Do not set it to zero. The value of rcoulomb (and likely anything related to
Coulombic interactions) is irrelevant because there will be no Coulombic
interactions if everything has zero charge.
-Justin
> Thank you kindly!
>
> Andrew DeYoung
> Carnegie Mellon University
>
> --
> Andrew DeYoung wrote:
>> Hi,
>>
>> I am running NVE simulations of a system of 1000 water molecules. I am
>> attempting to run analogous simulations on Gromacs 4.5.4 and also on
> another
>> MD program. Once finished, I will compare the results of the two
> programs.
>> It was suggested to me that to compare the results of the two programs in
> a
>> step-by-step manner, as a first step I should run the simulations with the
>> electrostatics (both short-range and long-range interactions) turned off.
>> In other words, run the simulations considering only Lennard-Jones
>> interactions. My question is, what is the best way to do this in Gromacs?
>>
>> A (very much more experienced) colleague suggested that I "turn off
>> electrostatics" by simply setting all the charges to zero. However, when
> I
>> look into the topology and structure files, I see some surprising (to me)
>> things -- the charges do not appear to be specified in the .top and .gro
>> files.
>>
>> To back up, I will first try to explain my system. I have a PDB file
>> specifying 1000 water molecules that someone gave me. This PDB file looks
>> like the following:
>>
>> TITLE water
>> REMARK THIS IS A SIMULATION BOX
>> CRYST1 47.100 47.100 47.100 90.00 90.00 90.00 P 1 1
>> MODEL 1
>> ATOM 1 OW SOL 1 -21.177 -21.195 -21.179 1.00 0.00
>> ATOM 2 HW1 SOL 1 -21.671 -21.195 -22.083 1.00 0.00
>> ATOM 3 HW2 SOL 1 -21.863 -21.195 -20.412 1.00 0.00
>> ATOM 4 OW SOL 2 -16.467 -21.195 -21.179 1.00 0.00
>> ATOM 5 HW1 SOL 2 -16.961 -21.195 -22.083 1.00 0.00
>> ATOM 6 HW2 SOL 2 -17.153 -21.195 -20.412 1.00 0.00
>> ...
>>
>> And so on. According to the PDB file format specification
>> (http://www.wwpdb.org/format/sect9.html), the last two columns in my file
>> correspond to the occupancies and temperature factors, respectively, of my
>> atom. So, evidently, partial charges on the oxygen and hydrogens are not
>> specified in this PDB file.
>>
>> Next, I execute the command:
>>
>> pdb2gmx -f water.pdb -p water.top -i water.itp -o water_processed.gro
>>
>> And, when prompted, I choose the OPLS-AA force field and the SPC/E water
>> model. Apparently, Gromacs recognizes from my PDB file that my system
>> contains water (and only water). Alternatively, I could have created a
>> topology file manually, without pdb2gmx, by simply creating a .top file
>> containing:
>> ------
>> #include "oplsaa.ff/forcefield.itp"
>> #include "oplsaa.ff/spce.itp"
>>
>> [ system ]
>> water
>>
>> [ molecules ]
>> SOL 1000
>> ------
>> and then simply feed grompp the topology file and the original PDB file
> when
>> preparing for a dynamics run (pdb2gmx, I understand, is meant for chains
>> rather than multiple molecules). But, in any case, I took the long road
> and
>> used pdb2gmx to create the .top topology and a .gro coordinate file.
>>
>> What surprises me (a novice in Gromacs and MD in general), though, is that
>> it seems that neither the .top file nor the .gro file contains any
> reference
>> to the charges on the oxygens and hydrogens in this water-only system.
> For
>> example, my .top file contains:
>> ------
>> ; Include forcefield parameters
>> #include "oplsaa.ff/forcefield.itp"
>>
>> ; Include water topology
>> #include "oplsaa.ff/spce.itp"
>>
>> #ifdef POSRES_WATER
>> ; Position restraint for each water oxygen
>> [ position_restraints ]
>> ; i funct fcx fcy fcz
>> 1 1 1000 1000 1000
>> #endif
>>
>> ; Include topology for ions
>> #include "oplsaa.ff/ions.itp"
>>
>> [ system ]
>> ; Name
>> water
>>
>> [ molecules ]
>> ; Compound #mols
>> SOL 1000
>> ------
>>
>> And, my .gro file contains, for example:
>> ------
>> water
>> 3000
>> 1HOH OW 1 -2.118 -2.119 -2.118
>> 1HOH HW1 2 -2.167 -2.119 -2.208
>> 1HOH HW2 3 -2.186 -2.119 -2.041
>> 2HOH OW 4 -1.647 -2.119 -2.118
>> 2HOH HW1 5 -1.696 -2.119 -2.208
>> 2HOH HW2 6 -1.715 -2.119 -2.041
>> ...
>> 4.71000 4.71000 4.71000
>> ------
>>
>> So, as far as I can tell, neither my .top nor my .gro file specify the
>> charges on the atoms. I guess that this means that the charges on the
> atoms
>> are only specified in the spce.itp file. When I find and open the
> spce.itp
>> file, I see:
>> ------
>> [ moleculetype ]
>> ; molname nrexcl
>> SOL 2
>>
>> [ atoms ]
>> ; nr type resnr residue atom cgnr charge mass
>> 1 opls_116 1 SOL OW 1 -0.8476
>> 2 opls_117 1 SOL HW1 1 0.4238
>> 3 opls_117 1 SOL HW2 1 0.4238
>>
>> #ifndef FLEXIBLE
>> [ settles ]
>> ; OW funct doh dhh
>> 1 1 0.1 0.16330
>>
>> [ exclusions ]
>> 1 2 3
>> 2 1 3
>> 3 1 2
>> #else
>> [ bonds ]
>> ; i j funct length force.c.
>> 1 2 1 0.1 345000 0.1 345000
>> 1 3 1 0.1 345000 0.1 345000
>>
>> [ angles ]
>> ; i j k funct angle force.c.
>> 2 1 3 1 109.47 383 109.47 383
>> #endif
>> ------
>>
>> So, there the charges are specified: -0.8476 on the oxygen and 0.4238 on
>> each hydrogen. Coming back to my original goal--to "turn off
>> electrostatics" by zeroing the charges--I could in principle just set the
>> charges in this .itp file to 0. However, on our cluster, this file
> appears
>> to be READ ONLY. So to change the charges to 0, I need to see if I can
> get
>> permission from our administrator.
>>
>
> That would be a terrible idea. It would affect anyone who uses that file,
> and
> then surprise! their entire simulation is worthless because the charges are
> not
> there :)
>
>> Based on my story above, my questions are:
>>
>> 1. Is there a better way to "turn off electrostatics"--without modifying
>> the spce.itp file?
>>
>
> There are several options:
>
> 1. Don't mess with the system-wide spce.itp file. Make your own local copy
> and
> alter it there. Adjust the call in the .top to some location where it will
> be
> found (which you can set with the "include" keyword in the .mdp file, if
> necessary).
>
> 2. Don't mess with topologies at all and instead use the free energy code
> with,
> e.g.:
>
> couple-moltype = SOL
> couple-lambda0 = vdw
> couple-moltype = (either vdw or none)
> init-lambda = 0
>
> Then don't both with any other lambda values, since you don't care about the
>
> free energy of some transformation of the system, just the energy value of
> the
> Hamiltonian under these conditions. This is a very roundabout way of
> accomplishing your task.
>
> 3. Use a simpler system like a noble gas, which has no charges to begin
> with.
>
>> 2. Why does it seem that charges are not specified in the .top and .gro
>> files? Is this just because I am boringly simulating only a solvent, and
>> the topology file tells Gromacs to look in spce.itp when simulating the
>> water molecules?
>>
>
> Charges *are* specified in the .top file, albeit indirectly. The #include
> statement is a C preprocessor term that literally means "cut and paste the
> contents of the included file right here." If you use grompp -pp to get the
>
> post-processed topology you will see. The #include statements are there for
>
> convenience, so you don't have to copy and paste a water topology for every
> single simulation, which would be horrendously annoying for most of us.
> Most
> other "normal" (i.e. protein-containing) topologies have all the information
> for
> the protein listed explicitly (see my tutorial, or probably any other, for a
>
> complete description).
>
> The .gro file does not contain charges because it is a coordinate file. It
> has
> atom and residue names and numbers, coordinates, and (optionally)
> velocities.
>
> -Justin
>
>> Thank you very much for your time!
>>
>> Andrew DeYoung
>> Carnegie Mellon University
>>
>
--
========================================
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
========================================
More information about the gromacs.org_gmx-users
mailing list