[gmx-users] Energy mimimization in AcOH box
Justin A. Lemkul
jalemkul at vt.edu
Sat Sep 12 15:36:10 CEST 2009
Aditi Borkar wrote:
> Thanx Justin.
> I had tried that too. Using both parameters I get the same outcome. I
> notice that increasing the emtol to 2000 (that effectively reduces the
> number of steps my system undergoes before converging) retains most of
> the structure. Should I change the dielectric constant of the solvent
> because my system is highly charged ( -744 e). The topology parameters
> from ASP side chain that I used are as follows:
I would take a step back and consider whether or not this model system is even
valid. A 1M solution of acetic acid would have a pH of roughly 2.4, meaning
that most of the molecules would be protonated; a 9M solution would probably be
at pH close to 0, meaning that having 744 negatively-charged molecules of
acetate probably does not come close to modeling reality.
If anything, you should probably be working with a net neutral system, since you
are probably experiencing a significant amount of charge repulsion in the
system. But strongly consider whether or not this model is appropriate.
> [ moleculetype ]
> ; Name nrexcl
> ACA 3
> [ atoms ]
> ; nr type resnr residue atom cgnr charge mass
> typeB chargeB massB
> 1 CH3 1 ACA CB 1 0 15.035 ; qtot 0
> 2 C 1 ACA CA 2 0.27 12.011
> ; qtot 0.27
> 3 OM 1 ACA OD1 2 -0.635 15.9994
> ; qtot -0.365
> 4 OM 1 ACA OD2 2 -0.635 15.9994 ; qtot -1
> [ bonds ]
> ; ai aj funct c0 c1 c2 c3
> 1 2 2 gb_26
> 2 3 2 gb_5
> 2 4 2 gb_5
> [ angles ]
> ; ai aj ak funct c0 c1 c2 c3
> 1 2 3 2 ga_21
> 1 2 4 2 ga_21
> 3 2 4 2 ga_37
> [ dihedrals ]
> ; ai aj ak al funct c0 c1
> c2 c3
> 1 4 3 2 2 gi_1
> -Aditi Borkar,
> Mumbai, India
> On Sat, Sep 12, 2009 at 5:35 PM, Justin A. Lemkul <jalemkul at vt.edu> wrote:
>> Aditi Borkar wrote:
>>> Dear All,
>>> I am simulating a protein in a box of 9M acetic acid. I have obtained
>>> the coordinates and gromacs topology for AcOH from PRODRG. The AcOH
>>> has net -1 charge and the topology file is as follows:
>> PRODRG topologies generally produce unsatisfactory charges. For acetic
>> acid, I would use the same charges found in an ASP side chain as a starting
>> point. You will note that your topology differs greatly from an ASP side
>> chain in any of the Gromos96 force fields.
>>> [ moleculetype ]
>>> ; Name nrexcl
>>> ACA 3
>>> [ atoms ]
>>> ; nr type resnr resid atom cgnr charge mass
>>> 1 CH3 1 ACA CAA 1 0.056 15.0350
>>> 2 C 1 ACA CAD 1 0.393 12.0110
>>> 3 OM 1 ACA OAC 1 -0.725 15.9994
>>> 4 OM 1 ACA OAB 1 -0.724 15.9994
>>> [ bonds ]
>>> ; ai aj fu c0, c1, ...
>>> 1 2 1 0.153 334720.0 0.153 334720.0 ; CAA CAD
>>> 2 3 1 0.125 418400.0 0.125 418400.0 ; CAD OAC
>>> 2 4 1 0.125 418400.0 0.125 418400.0 ; CAD OAB
>>> [ pairs ]
>>> ; ai aj fu c0, c1, ...
>>> [ angles ]
>>> ; ai aj ak fu c0, c1, ...
>>> 1 2 3 1 120.0 418.4 120.0 418.4 ; CAA CAD
>>> 1 2 4 1 120.0 418.4 120.0 418.4 ; CAA CAD
>>> 3 2 4 1 126.0 502.1 126.0 502.1 ; OAC CAD
>>> [ dihedrals ]
>>> ; ai aj ak al fu c0, c1, m, ...
>>> 2 1 4 3 2 0.0 1673.6 0.0 1673.6 ; imp CAD
>>> CAA OAB OAC
>>> When I had tried taking uncharged AcOH, during equilibration, the
>>> water and AcOH were completely segregating into two layers (like a
>>> biphasic system).
>>> My problem arises when I perform energy mimization on my solvated
>>> protein. With 300 steps of EM (emtol 1000 and emstep .01) the EM
>>> doesn't converge and more distrurbingly within these few steps, the
>>> secondary structure of my protein is completely lost. This problem
>>> does not occur when I do EM in box containing protein and pure water.
>>> The Fmax is displayed onto 2 TRP residues.
>> It is strange that this degree of change is occurring during EM. But I
>> wouldn't try anything else until you're sure you are starting with
>> reasonable parameters (see above).
>>> My em.mdp file is as follows:
>>> define = -DFLEXIBLE
>>> constraints = none
>>> integrator = steep
>>> dt = 0.002 ; ps !
>>> nsteps = 300
>>> nstlist = 10
>>> ns_type = grid
>>> rlist = 0.9
>>> coulombtype = PME
>>> rcoulomb = 0.9
>>> rvdw = 0.9
>>> fourierspacing = 0.12
>>> fourier_nx = 0
>>> fourier_ny = 0
>>> fourier_nz = 0
>>> pme_order = 4
>>> ewald_rtol = 1e-5
>>> optimize_fft = yes
>>> ; Energy minimizing stuff
>>> emtol = 1000.0
>>> emstep = 0.01
>>> Please help!
>>> Thank you!
>>> -Aditi Borkar,
>>> Tata Institute of Fundamental Research,
>>> gmx-users mailing list gmx-users at gromacs.org
>>> Please search the archive at http://www.gromacs.org/search before posting!
>>> Please don't post (un)subscribe requests to the list. Use the www
>>> interface or send it to gmx-users-request at gromacs.org.
>>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
>> Justin A. Lemkul
>> Ph.D. Candidate
>> ICTAS Doctoral Scholar
>> Department of Biochemistry
>> Virginia Tech
>> Blacksburg, VA
>> jalemkul[at]vt.edu | (540) 231-9080
>> gmx-users mailing list gmx-users at gromacs.org
>> Please search the archive at http://www.gromacs.org/search before posting!
>> Please don't post (un)subscribe requests to the list. Use the www interface
>> or send it to gmx-users-request at gromacs.org.
>> Can't post? Read http://www.gromacs.org/mailing_lists/users.php
Justin A. Lemkul
ICTAS Doctoral Scholar
Department of Biochemistry
jalemkul[at]vt.edu | (540) 231-9080
More information about the gromacs.org_gmx-users