[gmx-users] Protein in water/DMSO mixture with Amber

[Justin A. Lemkul]

Simone Cirri wrote:
> Hi all,
> to say the truth I already sent a message regarding this subject some
> months ago, but since so much time has passed I thought I'd open a new
> thread, also because the situation has quite evolved.
> 
> I'm trying to set up a simulation of a small (62 residues) protein in
> a mixture of water and DMSO with the AMBER99 forcefield. I found the
> parameters for such forcefield in the literature and used them to
> build a dmso.itp. Following yuor previous suggestions, I also obtained
> a dmso.gro file downloading a .sdf file from Pubchem and converting it
> to a .pdb and eventually to a .gro.
> Then I followed the mini-tutorial in the Gromacs site on simulations
> of mixed solvents: I added some DMSO molecules with genbox -ci -nmol
> and then filled the box with TIP3P water molecules. I also added 9 Cl-
> ions and, judging from what I could see with a visualization program,
> everything was fine.
> The problem came when I tried to run the energy minimization. I got the error:
> 
> Stepsize too small, or no change in energy.
> Converged to machine precision,
> but not to the requested precision Fmax < 100
> 
> Double precision normally gives you higher accuracy.
> You might need to increase your constraint accuracy, or turn
> off constraints alltogether (set constraints = none in mdp file)
> 
> after only 5 minimization steps. I searched in the GROMACS site and in
> this mailing list as well to find out what the error meant and I read
> that this error can be sometimes ignored if the potential energy of
> the system is already low enough to run the simulation... but my Epot
> is 7.6367096e+16 . To be completely sure, I tried to run an MD
> simulation, but it immediately crashed with a "water molecule cannot
> be settled" error, which does not surprise me.
> 
> What did I do wrong? Do I have to assume that there is something wrong
> in my dmso.itp file or in the way I ran the minimization?
> Thank you very much; I'm pasting my dmso.itp and em.mdp here.
> 

You could have some unresolvable atomic overlap in the system, so visual 
inspection is your friend there.  You can also try using a larger value of 
emstep, something like 0.01, otherwise you could just be stuck in a high-energy 
region from which your structure won't escape using the tiny emstep you have set.

-Justin

> ###DMSO.ITP###
> 
> [ moleculetype ]
> ; molname  nrexcl
> DMSO	2
> 
> [ atoms ]
> #ifdef _FF_AMBER99
> ;   nr    type   	resnr  residue atom    cgnr    charge  	mass
>     1     amber99_47	1	DMSO	S	1	0.3155  32.06000                 ;
> amber S  type
>     2     amber99_41	1	DMSO	O	1	-0.5205 16.00000                 ;
> amber O  type
>     3	  amber99_11	1	DMSO	CT1	1	-0.3244 12.01000                 ;
> amber C  type
>     4     amber99_19    1	DMSO	HT11	1	0.1423  1.00800
> ; amber H  type
>     5     amber99_19    1       DMSO    HT12      1       0.1423
> 1.00800                 ; amber H  type
>     6     amber99_19    1       DMSO    HT13      1       0.1423
> 1.00800                 ; amber H  type
>     7	  amber99_11    1       DMSO    CT2      1       -0.3244
> 12.01000                 ; amber C  type
>     8     amber99_19    1	DMSO	HT21	1	0.1423  1.00800
> ; amber H  type
>     9     amber99_19    1       DMSO    HT22      1       0.1423
> 1.00800                 ; amber H  type
>     10    amber99_19    1       DMSO    HT23      1       0.1423
> 1.00800                 ; amber H  type
> #endif
> 
> [ bonds ]
> ;  ai    aj funct          b0       kb
>     1     2     1 	0.153	    238602
>     1     3     1 	0.181	    95022
>     3     4     1 	0.109	    142324
>     3	  5     1	0.109       142324
>     3     6     1	0.109       142324
>     1     7     1	0.181       95022
>     7     8     1	0.109       142324
>     7     9     1	0.109       142324
>     7     10    1	0.109       142324
> 
> [ angles ]
> ;  ai    aj    ak   funct        theta	     kt
> 3	1	2	1	 106.75      334.88
> 3	1	7	1	 97.40	     259.53
> 2	1	7	1	 106.75	     334.88
> 4       3       1       1	 109.50	     209.30
> 6	3	1	1	 109.50      209.30
> 5	3	1	1	 109.50      209.30
> 5	3	4	1	 109.50	     146.51
> 5	3	6	1	 109.50	     146.51
> 4	3	6	1	 109.50	     146.51
> 9	7	10	1	 109.50	     146.51
> 8	7	9	1	 109.50	     146.51
> 8	7	10	1	 109.50	     146.51
> 9	7	1	1	 109.50      209.30
> 10	7	1	1	 109.50      209.30
> 8	7	1	1	 109.50      209.30
> 
> [ dihedrals ]
> ;  ai  aj      ak      al   funct     xi         kxi	pn
> 2	1	3	6      1      180        3.558	2
> 2	1	3	4      1      180        3.558	2
> 2	1	3	5      1      180        3.558	2
> 2	1	7	8      1      180        3.558	2
> 2	1	7	9      1      180        3.558	2
> 2	1	7	10     1      180        3.558	2
> 3	1	7	8      1      180        3.558	2
> 3	1	7	9      1      180        3.558	2
> 3	1	7	10     1      180        3.558	2
> 7	1	3	4      1      180        3.558	2
> 7	1	3	5      1      180        3.558	2
> 7	1	3	6      1      180        3.558	2
> 
> 
> ###EM.MDP###
> 
> title               =  1HAA
> cpp                 =  /lib/cpp
> define              =  -DFLEX_SPC
> constraints         =  none
> integrator          =  steep
> dt                  =  0.002    ; ps !
> nsteps              =  10000
> 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		=  6
> ewald_rtol		=  1e-5
> optimize_fft		=  yes
> ;
> ;       Energy minimizing stuff
> ;
> emtol               =  100.0
> emstep              =  0.00001

-- 
========================================

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

========================================