[gmx-users] Crystallization of NaCl saturated aqueous solution under high temperature of low pressure.

[Justin Lemkul]

On 9/11/17 1:53 AM, Yanke Peng wrote:
> Hi, Justin, Thanks very much for your advice.
> I think I should try a low concentration NaCl to verify if the system has
> the chance to be energy minimizated.
>
> Here are the exact commands that I used.
> 1. solution definition
> #water box and topology file generation
> gmx solvate -box 2.48 2.48 2.56 -o water.gro -maxsol 352

Does this actually fill the box appropriately, and do you actually get 
all 352 waters?

> gmx pdb2gmx -f water.gro -p water.top
> #get water.tpr
> gmx grompp -f water.mdp -c water.gro -p water.top -o water.tpr
>    (#content of water.mdp file: integrator=steep)
> #get solution.gro
> gmx genion -s water.tpr -np 32 -pname NA -nn 32 -nname CL -o solution.gro
> -p water.top
>   (after 32Na+ and 5Cl+ was replaced. The fatal error shows "No more
> replaceable solvent")

This should be your first indication of a problem.  If you actually have 
352 waters, there should be no issue adding 64 ions.  Visualize the 
structures at each step, and verify the contents of the topology.  If 
you ask for 32 of each ion, you should get it.  If you don't, stop and 
re-evaluate what you have done.

-Justin

> 2. energy minimizaiton
> # em.mdp
> (integrator = steep
> nsteps     = 10000
>
> cutoff-scheme = Verlet
> ns_type       = grid
> pbc           = xyz
>
> coulombtype      = PME
> coulomb-modifier = Potential-shift-Verlet
> rcoulomb         = 1.2
>
> vdwtype      = Cut-off
> vdw-modifier = Potential-shift-Verlet
> rvdw         = 1.2
> DispCorr     = no
>
> constraints = hbonds)
>
> #get em.tpr and do the EM
> gmx grompp -f em.mdp -c solution.gro -p solution.top -o em.tpr
> gmx mdrun -v -deffnm em
> //after this command, the information shows as below.
> Steepest Descents:
>     Tolerance (Fmax)   =  1.00000e+01
>     Number of steps    =        10000
> Step=    0, Dmax= 1.0e-02 nm, Epot= -5.99726e+03 Fmax= 3.44081e+04, atom=
> 917
> Step=    1, Dmax= 1.0e-02 nm, Epot= -7.14350e+03 Fmax= 1.48023e+04, atom=
> 907
> Step=    2, Dmax= 1.2e-02 nm, Epot= -8.56126e+03 Fmax= 7.60131e+03, atom=
> 907
> ....(step 3 to 11 was abridged here)
> Step=   12, Dmax= 7.4e-02 nm, Epot= -2.39035e+04 Fmax= 4.10568e+03, atom=
> 878
> Step=   13, Dmax= 8.9e-02 nm, Epot= -2.46762e+04 Fmax= 5.66839e+03, atom=
> 878
>
> Back Off! I just backed up step14b.pdb to ./#step14b.pdb.2#
>
> fatal error shows "step 14: Water molecule starting at atom 22 can not be
> settled."
>
> # *then, I changed the integrator=md. No error information showed up. but
> the em.trr file is empty.*
>
> 3. equilibration
> #contents of eq-nvt.mdp
> (integrator = md
> dt         = 0.001  ;1 fs
> nsteps     = 10000  ;10 ps
>
> cutoff-scheme = Verlet
> ns_type       = grid
> pbc           = xyz
>
> coulombtype      = PME
> coulomb-modifier = Potential-shift-Verlet
> rcoulomb         = 1.2
>
> vdwtype      = Cut-off
> vdw-modifier = Potential-shift-Verlet
> rvdw         = 1.2
> DispCorr     = no
>
> constraints = hbonds
>
> tcoupl  = v-rescale
> tau-t   = 1.0
> ref-t   = 298
> tc-grps = System
>
> gen_vel  = yes
> gen_temp = 298)
>
> # do the equilibration
> Command line:
>    gmx grompp -f eq-nvt.mdp -c em.tpr -t *em.trr* -p solution.top -o
> eq-nvt.tpr
>
>
> Back Off! I just backed up mdout.mdp to ./#mdout.mdp.2#
> Setting the LD random seed to 3348907021 <(334)%20890-7021>
> Generated 2211 of the 2211 non-bonded parameter combinations
> Generating 1-4 interactions: fudge = 0.5
> Generated 2211 of the 2211 1-4 parameter combinations
> Excluding 2 bonded neighbours molecule type 'SOL'
> turning H bonds into constraints...
> Excluding 1 bonded neighbours molecule type 'NA'
> turning H bonds into constraints...
> Excluding 1 bonded neighbours molecule type 'CL'
> turning H bonds into constraints...
> Reading file ../1-em/em.tpr, VERSION 5.1.2 (single precision)
> Reading file ../1-em/em.tpr, VERSION 5.1.2 (single precision)
> Setting gen_seed to 345062238
> Velocities were taken from a Maxwell distribution at 298 K
> Removing all charge groups because cutoff-scheme=Verlet
> Analysing residue names:
> There are:   288      Water residues
> There are:    64        Ion residues
> Analysing residues not classified as Protein/DNA/RNA/Water and splitting
> into groups...
> Number of degrees of freedom in T-Coupling group System is 1917.00
> Determining Verlet buffer for a tolerance of 0.005 kJ/mol/ps at 298 K
> Calculated rlist for 1x1 atom pair-list as 1.202 nm, buffer size 0.002 nm
> Set rlist, assuming 4x4 atom pair-list, to 1.200 nm, buffer size 0.000 nm
> Note that mdrun will redetermine rlist based on the actual pair-list setup
> Reading Coordinates and Box size from old trajectory
> Will read whole trajectory
> Velocities generated: ignoring velocities in input trajectory
> Last frame         -1 time    0.000
>
> -------------------------------------------------------
> Program gmx grompp, VERSION 5.1.2
> Source code file: /build/gromacs-z6bPBg/gromacs-
> 5.1.2/src/gromacs/gmxpreprocess/grompp.c, line: 778
>
> Fatal error:
> *Number of atoms in Topology is not the same as in Trajectory*
> (the input file em.trr is empty as shown before, so the failure of this
> step is predictable)
>
>
> On Sun, Sep 10, 2017 at 5:28 PM, Yanke Peng <pengx288 at umn.edu> wrote:
>
>> Hello to all,
>>
>> I am trying to simulate the crystallization process of NaCl saturated
>> aqueous solution by the evaporation of solvent - water molecules under high
>> temperature or low pressure.
>>
>> Basic settings: 352 water molecules in a 2.5*2.5*2.5 nm box;  For a
>> saturated solution, 32 Na+ and 32 Cl- should be introduced by replacing 64
>> water molecules.
>>
>> Below are the problems I encountered.
>>
>> 1.     gmx genion only could replace 37 water molecules with 32 Na+ and 5
>> Cl-.  Does this means gmx genion couldn’t generate saturated solution
>> sometimes? So, I used Avogadro to generate a coordinate file of saturated
>> NaCl solution and then generated a .gro file by the coordinate file.
>>
>> 2.     During the energy minimization,  I chose integrator=steep in em.mdp
>> and after 13 steps, the process stopped.  Here is the error information.
>>
>> Fatal error:
>>
>> step 14: Water molecule starting at atom 22 can not be settled.
>> Check for bad contacts and/or reduce the timestep if appropriate.
>>
>> 3. I changed the integrator=md and no error information showed up during
>> the -mdrun.  However, the em.trr file is empty.  So when I rum the grompp
>> during the equilibration NVT, the error information said "Number of atoms
>> in Topology is not the same as in Trajectory".
>>
>> I am trying to understand the mechanism of gromacs by learning the manual
>> thoroughly.
>>
>> Could you give me a clue to solve these problems? I really appreciate any
>> help from you!
>>
>>
>> <pengx288 at umn.edu>
>>
>>
>>

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

Justin A. Lemkul, Ph.D.
Assistant Professor
Virginia Tech Department of Biochemistry

303 Engel Hall
340 West Campus Dr.
Blacksburg, VA 24061

jalemkul at vt.edu | (540) 231-3129
http://www.biochem.vt.edu/people/faculty/JustinLemkul.html

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