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

[Yanke Peng]

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
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")
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>
>
>
>